diff --git a/.vscode/launch.json b/.vscode/launch.json index 61dfa9517..6340699ff 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -9,7 +9,7 @@ "console": "integratedTerminal", "args": [ "-p", - "${workspaceFolder}/sharc/input/parameters.yaml" + "${workspaceFolder}/sharc/campaigns/ast_mss_d2d_to_imt_cpe/input/parameter_ast_mss_d2d_to_imt_cpe_30exclusion_45beam_elev_0.0power_backoff_20load_factor_P619prop_imt-ue_imt.upto-1GHz.single-bs.urban-macro-bs_system-4.698-960MHz-block2.690km-antenna-update.yaml" ] } ] diff --git a/sharc/antenna/antenna_f1245_fs.py b/sharc/antenna/antenna_f1245_fs.py new file mode 100644 index 000000000..de554a19d --- /dev/null +++ b/sharc/antenna/antenna_f1245_fs.py @@ -0,0 +1,196 @@ + +# -*- coding: utf-8 -*- +""" +Created on Wed Apr 4 17:08:00 2018 +@author: Calil +""" +import matplotlib.pyplot as plt +from sharc.antenna.antenna import Antenna +from sharc.parameters.imt.parameters_imt import ParametersImt +from sharc.parameters.imt.parameters_imt import ParametersImt +from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt +from sharc.parameters.parameters_antenna import ParametersAntenna +import numpy as np +import math + + +class Antenna_f1245_fs(Antenna): + """Class that implements the ITU-R F.1245 antenna pattern for fixed + satellite service earth stations.""" + + def __init__(self, param: ParametersImt): + super().__init__() + self.peak_gain = param.gain + lmbda = 3e8 / (param.frequency * 1e6) + self.d_lmbda = param.diameter / lmbda + self.g_l = 2 + 15 * math.log10(self.d_lmbda) + self.phi_m = (20 / self.d_lmbda) * math.sqrt(self.peak_gain - self.g_l) + self.phi_r = 12.02 * math.pow(self.d_lmbda, -0.6) + + def calculate_gain(self, *args, **kwargs) -> np.array: + """ + Calculate the antenna gain for the given parameters. + + Parameters + ---------- + *args : tuple + Positional arguments (not used). + **kwargs : dict + Keyword arguments, expects 'phi_vec', 'theta_vec', and 'beams_l'. + + Returns + ------- + np.array + Calculated antenna gain values. + """ + # phi_vec = np.absolute(kwargs["phi_vec"]) + # theta_vec = np.absolute(kwargs["theta_vec"]) + # beams_l = np.absolute(kwargs["beams_l"]) + off_axis = np.absolute(kwargs["off_axis_angle_vec"]) + if self.d_lmbda > 100: + gain = self.calculate_gain_greater(off_axis) + else: + gain = self.calculate_gain_less(off_axis) + # idx_max_gain = np.where(beams_l == -1)[0] + # gain = self.peak_gain + return gain + + def calculate_gain_greater(self, phi: float) -> np.array: + """ + For frequencies in the range 1 GHz to about 70 GHz, in cases where the + ratio between the antenna diameter and the wavelength is GREATER than + 100, this method should be used. + Parameter + --------- + phi : off-axis angle [deg] + Returns + ------- + a numpy array containing the gains in the given angles + """ + gain = np.zeros(phi.shape) + idx_0 = np.where(phi < self.phi_m)[0] + gain[idx_0] = self.peak_gain - 2.5e-3 * \ + np.power(self.d_lmbda * phi[idx_0], 2) + phi_thresh = max(self.phi_m, self.phi_r) + idx_1 = np.where((self.phi_m <= phi) & (phi < phi_thresh))[0] + gain[idx_1] = self.g_l + idx_2 = np.where((phi_thresh <= phi) & (phi < 48))[0] + gain[idx_2] = 29 - 25 * np.log10(phi[idx_2]) + idx_3 = np.where((48 <= phi) & (phi <= 180))[0] + gain[idx_3] = -13 + return gain + + def calculate_gain_less(self, phi: float) -> np.array: + """ + For frequencies in the range 1 GHz to about 70 GHz, in cases where the + ratio between the antenna diameter and the wavelength is LESS than + or equal to 100, this method should be used. + Parameter + --------- + phi : off-axis angle [deg] + Returns + ------- + a numpy array containing the gains in the given angles + """ + gain = np.zeros(phi.shape) + idx_0 = np.where(phi < self.phi_m)[0] + gain[idx_0] = self.peak_gain - 0.0025 * \ + np.power(self.d_lmbda * phi[idx_0], 2) + idx_1 = np.where((self.phi_m <= phi) & (phi < 48))[0] + gain[idx_1] = 39 - 5 * \ + math.log10(self.d_lmbda) - 25 * np.log10(phi[idx_1]) + idx_2 = np.where((48 <= phi) & (phi < 180))[0] + gain[idx_2] = -3 - 5 * math.log10(self.d_lmbda) + return gain + + def add_beam(self, phi: float, theta: float): + """ + Add a new beam to the antenna. + + Parameters + ---------- + phi : float + Azimuth angle in degrees. + theta : float + Elevation angle in degrees. + """ + self.beams_list.append((phi, theta)) + + def calculate_off_axis_angle(self, Az, b): + """ + Calculate the off-axis angle between the main beam and a given direction. + + Parameters + ---------- + Az : float or np.array + Azimuth angle(s) in degrees. + b : float or np.array + Elevation angle(s) in degrees. + + Returns + ------- + float or np.array + Off-axis angle(s) in degrees. + """ + Az0 = self.beams_list[0][0] + a = 90 - self.beams_list[0][1] + C = Az0 - Az + off_axis_rad = np.arccos( + np.cos( + np.radians(a)) * + np.cos( + np.radians(b)) + + np.sin( + np.radians(a)) * + np.sin( + np.radians(b)) * + np.cos( + np.radians(C)), + ) + off_axis_deg = np.degrees(off_axis_rad) + return off_axis_deg + + +if __name__ == '__main__': + off_axis_angle_vec = np.linspace(0.1, 180, num=1001) + # initialize antenna parameters + param = ParametersAntenna() + param.frequency = 2155 + param_gt = ParametersAntennaImt() + param.gain = 33.1 + param.diameter = 2 + antenna_gt = Antenna_f1245_fs(param) + antenna_gt.add_beam(0, 0) + gain_gt = antenna_gt.calculate_gain( + off_axis_angle_vec=off_axis_angle_vec, + ) + param.diameter = 3 + antenna_gt = Antenna_f1245_fs(param) + gain_gt_3 = antenna_gt.calculate_gain( + off_axis_angle_vec=off_axis_angle_vec, + ) + param.diameter = 1.8 + antenna_gt = Antenna_f1245_fs(param) + gain_gt_18 = antenna_gt.calculate_gain( + off_axis_angle_vec=off_axis_angle_vec, + ) + + fig = plt.figure( + figsize=(8, 7), facecolor='w', + edgecolor='k', + ) # create a figure object + plt.semilogx(off_axis_angle_vec, gain_gt, "-b", label="$f = 10.7$ $GHz,$ $D = 2$ $m$") + plt.semilogx(off_axis_angle_vec, gain_gt_3, "-y", label="$f = 10.7$ $GHz,$ $D = 3$ $m$") + plt.semilogx(off_axis_angle_vec, gain_gt_18, "-g", label="$f = 10.7$ $GHz,$ $D = 1.8$ $m$") + + plt.title("ITU-R F.1245 antenna radiation pattern") + plt.xlabel(r"Off-axis angle $\phi$ [deg]") + plt.ylabel("Gain relative to $G_m$ [dB]") + plt.legend(loc="lower left") + # plt.xlim((phi[0], phi[-1])) + plt.ylim((-20, 50)) + # ax = plt.gca() + # ax.set_yticks([-30, -20, -10, 0]) + # ax.set_xticks(np.linspace(1, 9, 9).tolist() + np.linspace(10, 100, 10).tolist()) + plt.grid() + plt.show() diff --git a/sharc/antenna/antenna_factory.py b/sharc/antenna/antenna_factory.py index 051e1f332..e6e33c6ae 100644 --- a/sharc/antenna/antenna_factory.py +++ b/sharc/antenna/antenna_factory.py @@ -12,6 +12,7 @@ from sharc.antenna.antenna_s580 import AntennaS580 from sharc.antenna.antenna_s1528 import AntennaS1528 from sharc.antenna.antenna_s1855 import AntennaS1855 +from sharc.antenna.antenna_f1245_fs import Antenna_f1245_fs from sharc.antenna.antenna_s1528 import AntennaS1528, AntennaS1528Leo, AntennaS1528Taylor from sharc.antenna.antenna_beamforming_imt import AntennaBeamformingImt @@ -49,6 +50,8 @@ def create_antenna( return AntennaS465(antenna_params.itu_r_s_465_modified) case "ITU-R S.1855": return AntennaS1855(antenna_params.itu_r_s_1855) + case "ITU-R F.1245_fs": + return Antenna_f1245_fs(antenna_params.itu_r_f_1245_fs) case "ITU-R Reg. RR. Appendice 7 Annex 3": return AntennaReg_RR_A7_3(antenna_params.itu_reg_rr_a7_3) case "MSS Adjacent": @@ -60,6 +63,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 { @@ -88,6 +97,14 @@ def create_n_antennas( antennas[i] = AntennaFactory.create_antenna( antenna_params, azimuth[i], elevation[i], ) + elif antenna_params.pattern == "Antenna System 4": + antenna_pattern_high = AntennaS1528(antenna_params.antenna_system_4.antenna_parameters_high) + antenna_pattern_low = AntennaS1528(antenna_params.antenna_system_4.antenna_parameters_low) + for i in range(n_stations): + if antenna_params.antenna_system_4.beam_ground_elev_angles[i] >= 50: + antennas[i] = antenna_pattern_high + else: + antennas[i] = antenna_pattern_low else: # some antennas don't need azimuth and elevation at all # this makes it much faster diff --git a/sharc/parameters/antenna/parameters_antenna_system4.py b/sharc/parameters/antenna/parameters_antenna_system4.py new file mode 100644 index 000000000..340bbf65f --- /dev/null +++ b/sharc/parameters/antenna/parameters_antenna_system4.py @@ -0,0 +1,52 @@ +# 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 + ) + ) + + # Elevation angles (degrees) separating high and low elevation beams + # This is set during simulation runtime. + beam_ground_elev_angles: list[float] = None + + 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_imt_mss_dc.py b/sharc/parameters/imt/parameters_imt_mss_dc.py index c14079f4f..feffbb32c 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): + self.zones[i].power_backoff_db = 0.0 # Default to 0 dB if not set - we make sure it's not applied later. + + @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 a8fc8a0e0..946d247bb 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,9 @@ 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", + "ITU-R F.1245_fs"] # chosen antenna radiation pattern pattern: typing.Literal["OMNI", @@ -44,7 +47,9 @@ 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", + "ITU-R F.1245_fs"] = None # antenna gain [dBi] gain: float = None @@ -94,6 +99,25 @@ class ParametersAntenna(ParametersBase): default_factory=ParametersAntennaS672, ) + antenna_system_4: ParametersAntennaSystem4 = field( + default_factory=ParametersAntennaSystem4, + ) + + @dataclass + class ParametersAntennaRF1245(ParametersBase): + gain: float = -25 + diameter: float = None + frequency: float = None + + def validate(self, ctx): + """ do a sanity check of the parameters """ + if None in [self.gain, self.diameter, self.frequency]: + raise ValueError(f"{ctx}.antenna_3_dB should be set to a number") + + itu_r_f_1245_fs: ParametersAntennaRF1245 = field( + default_factory=ParametersAntennaRF1245, + ) + def set_external_parameters(self, **kwargs): """ Set external parameters for all sub-parameters of the antenna. @@ -198,6 +222,10 @@ 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 "ITU-R F.1245_fs": + self.itu_r_f_1245_fs.validate(f"{ctx}.itu_r_f_1245_fs") 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 c566269c0..7b0e5961c 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/post_processor.py b/sharc/post_processor.py index 9a8f67431..59b4d4a5c 100644 --- a/sharc/post_processor.py +++ b/sharc/post_processor.py @@ -159,6 +159,10 @@ class PostProcessor: "title": "[IMT] UL throughput", "x_label": "Throughput [bits/s/Hz]", }, + "imt_ul_interf_power": { + "title": "[IMT] UL interference power", + "x_label": "Interference Power [dB]", + }, "imt_path_loss": { "title": "[IMT] path loss", "x_label": "Path loss [dB]", @@ -223,6 +227,10 @@ class PostProcessor: "title": "[IMT] DL throughput", "x_label": "Throughput [bits/s/Hz]", }, + "imt_dl_interf_power": { + "title": "[IMT] DL interference power", + "x_label": "Interference Power [dB]", + }, "system_ul_interf_power": { "title": "[SYS] system interference power from IMT UL", "x_label": "Interference Power [dBm/BMHz]", @@ -247,12 +255,12 @@ class PostProcessor: "x_label": "Transmit power [dBm]", "title": "[IMT] DL transmit power", }, - "imt_dl_pfd_external": { - "title": "[IMT] DL external Power Flux Density (PFD) ", + "imt_dl_pfd": { + "title": "[IMT] DL Power Flux Density (PFD) ", "x_label": "PFD [dBW/m²/MHz]", }, - "imt_dl_pfd_external_aggregated": { - "title": "[IMT] Aggregated DL external Power Flux Density (PFD)", + "imt_dl_aggregated_pfd": { + "title": "[IMT] DL Aggregated Power Flux Density (PFD) ", "x_label": "PFD [dBW/m²/MHz]", }, # these ones were not plotted already, so will continue to not be diff --git a/sharc/results.py b/sharc/results.py index 92bd005da..926cfe73b 100644 --- a/sharc/results.py +++ b/sharc/results.py @@ -42,6 +42,8 @@ def __init__(self): self.imt_ul_tput_ext = SampleList() # Throughput [bits/s/Hz] self.imt_ul_tput = SampleList() + # System to IMT interference power [dB] - UL + self.imt_ul_interf_power = SampleList() self.imt_path_loss = SampleList() self.imt_coupling_loss = SampleList() @@ -81,10 +83,13 @@ def __init__(self): self.imt_dl_tput_ext = SampleList() # Throughput [bits/s/Hz] self.imt_dl_tput = SampleList() + # System to IMT interference power [dB] - DL + self.imt_dl_interf_power = SampleList() - # PFD - self.imt_dl_pfd_external = SampleList() - self.imt_dl_pfd_external_aggregated = SampleList() + # PFD at the IMT receiver [dBW/m²/MHz] + self.imt_dl_pfd = SampleList() + self.imt_dl_aggregated_pfd = SampleList() + self.imt_ul_pfd = SampleList() self.system_ul_coupling_loss = SampleList() self.system_ul_interf_power = SampleList() diff --git a/sharc/simulation.py b/sharc/simulation.py index 07e44a800..a7d24adf3 100644 --- a/sharc/simulation.py +++ b/sharc/simulation.py @@ -105,6 +105,7 @@ def __init__(self, parameters: Parameters, parameter_file: str): self.imt_bs_antenna_gain = list() self.imt_ue_antenna_gain = list() self.system_imt_antenna_gain = list() + self.system_imt_pfd = list() self.imt_system_antenna_gain = list() self.imt_system_antenna_gain_adjacent = list() self.imt_system_path_loss = list() @@ -745,6 +746,31 @@ def calculate_bw_weights( return overlap + def calculate_system_to_imt_pfd(self, imt_station: StationManager): + """Calculates the Power Flux Density (PFD) at the IMT station + + PFD formula (dBW/m²/MHz) + PFD = EIRP - 10log10(4π) - 20log10(distance) + + Parameters + ---------- + imt_station : StationManager + A StationManager object with IMT stations + """ + # EIRP in dBW/MHz + if imt_station.station_type is StationType.IMT_BS: + sys_to_imt_antenna_gain = \ + self.system_imt_antenna_gain[:, np.arange(0, + self.system_imt_antenna_gain.shape[1] - 1, + self.parameters.imt.ue.k)] + else: + sys_to_imt_antenna_gain = self.system_imt_antenna_gain + + eirp_dBW_MHz = self.system.tx_power_density[:, np.newaxis] + 60 + \ + sys_to_imt_antenna_gain + dist_sys_to_imt = self.system.geom.get_3d_distance_to(imt_station.geom) + self.system_imt_pfd = eirp_dBW_MHz - 10.992098640220963 - 20 * np.log10(dist_sys_to_imt) + def plot_scenario(self): """ Plot the current simulation scenario, including network topology and user equipments. diff --git a/sharc/simulation_downlink.py b/sharc/simulation_downlink.py index 58e9231be..0ebb667fe 100644 --- a/sharc/simulation_downlink.py +++ b/sharc/simulation_downlink.py @@ -12,7 +12,7 @@ from sharc.simulation import Simulation from sharc.parameters.parameters import Parameters from sharc.station_factory import StationFactory -from sharc.parameters.constants import BOLTZMANN_CONSTANT +from sharc.parameters.constants import BOLTZMANN_CONSTANT, SPEED_OF_LIGHT from sharc.propagation.propagation_path import PropagationPath warn = warnings.warn @@ -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 * \ @@ -408,32 +409,8 @@ def calculate_sinr_ext(self): self.ue.inr[ue] = self.ue.ext_interference[ue] - \ self.ue.thermal_noise[ue] - # Calculate PFD at the UE - - # Distance from each system transmitter to each UE receiver (in meters) - dist_sys_to_imt = self.system.geom.get_3d_distance_to( - self.ue.geom) # shape: [n_tx, n_ue] - - # EIRP in dBW/MHz per transmitter - eirp_dBW_MHz = self.param_system.tx_power_density + \ - 60 + self.system_imt_antenna_gain - - # PFD formula (dBW/m²/MHz) - # PFD = EIRP - 10log10(4π) - 20log10(distance) - # Store the PFD for each transmitter and each UE - self.ue.pfd_external = eirp_dBW_MHz - \ - 10.992098640220963 - 20 * np.log10(dist_sys_to_imt) - - # Total PFD per UE (sum of PFDs from each transmitter) - # Convert PFD from dB to linear scale (W/m²/MHz) - pfd_linear = 10 ** (self.ue.pfd_external / 10) - # Sum PFDs from all transmitters for each UE (axis=0 assumes shape - # [n_tx, n_ue]) - sys_active = np.where(self.system.active)[0] - # FIXME: consider only correct paths here - pfd_agg_linear = np.sum(pfd_linear[sys_active], axis=0) - # Convert back to dBW - self.ue.pfd_external_aggregated = 10 * np.log10(pfd_agg_linear) + # Calculate PFD at the UE + self.calculate_system_to_imt_pfd(self.ue) def calculate_external_interference(self): """ @@ -472,14 +449,6 @@ def calculate_external_interference(self): self.param_system.frequency, ) - interference = self.bs.tx_power[frst_bs] - pow_coch = 10 * np.log10( - weights * np.power( - 10, - 0.1 * interference, - ), - ) - if self.adjacent_channel: # Calculate how much power is emitted in the adjacent channel: if self.parameters.imt.adjacent_ch_emissions == "SPECTRAL_MASK": @@ -564,6 +533,12 @@ def calculate_external_interference(self): ) ] if self.co_channel: + pow_coch = 10 * np.log10( + weights * np.power( + 10, + 0.1 * self.bs.tx_power[bs], + ), + ) rx_interference += np.sum( 10 ** (0.1 * (pow_coch - self.coupling_loss_imt_system[active_beams, system_interfering])) ) @@ -652,8 +627,8 @@ def collect_results(self, write_to_file: bool, snapshot_number: int): ) self.add_system_imt_interaction_attr_to_results( "DL", - self.ue.pfd_external, - "imt_dl_pfd_external", + self.system_imt_pfd, + "imt_dl_pfd", ) self.add_system_imt_interaction_attr_to_results("DL", "system_imt_antenna_gain") @@ -673,6 +648,29 @@ def collect_results(self, write_to_file: bool, snapshot_number: int): for bs in bs_active: ue = self.link[bs] + ############################################# + # Experimental PFD calculation at UE based on received power + if self.parameters.imt.ue.antenna.pattern == "ARRAY": + Gr = self.parameters.imt.ue.antenna.array.element_max_g + if self.parameters.imt.ue.antenna.array.element_pattern != "FIXED": + raise NotImplementedError("PFD with non-FIXED element pattern is not implemented") + elif self.parameters.imt.ue.antenna.pattern != "OMNI": + raise NotImplementedError("This antenna pattern is not implemented for PFD calc") + else: + Gr = self.parameters.imt.ue.antenna.gain + + # TODO: check if also remove polarization_loss + L = self.parameters.imt.ue.ohmic_loss \ + + self.parameters.imt.ue.body_loss + wavelen = SPEED_OF_LIGHT / (self.parameters.imt.frequency * 1e6) + # self.ue.ext_interference[ue] is already without noise + # and after coupling loss + # WARNING: overwriting stuff + # self.ue.ext_interference is in dBm + self.ue.pfd_external_aggregated[ue] = self.ue.ext_interference[ue] - 30 - Gr + L - \ + 10 * np.log10(wavelen**2 / (4 * np.pi)) - 10 * np.log10(self.ue.bandwidth[ue]) + ############################################# + if not self.parameters.imt.imt_dl_intra_sinr_calculation_disabled: self.results.imt_path_loss.extend(self.path_loss_imt[bs, ue]) self.results.imt_coupling_loss.extend( @@ -708,8 +706,9 @@ def collect_results(self, write_to_file: bool, snapshot_number: int): ) self.results.imt_dl_inr.extend(self.ue.inr[ue].tolist()) - self.results.imt_dl_pfd_external_aggregated.extend( - self.ue.pfd_external_aggregated[ue].tolist()) + self.results.imt_dl_interf_power.extend( + (self.ue.ext_interference[ue] - 30).tolist(), + ) self.results.imt_dl_tx_power.extend(self.bs.tx_power[bs].tolist()) diff --git a/sharc/simulation_uplink.py b/sharc/simulation_uplink.py index ebb307189..c50356e42 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)) @@ -349,6 +348,9 @@ def calculate_sinr_ext(self): self.bs.inr[bs] = self.bs.ext_interference[bs] - \ self.bs.thermal_noise[bs] + # Calculate PFD at the BS + self.calculate_system_to_imt_pfd(self.bs) + def calculate_external_interference(self): """ Calculates interference that IMT system generates on other system @@ -581,6 +583,12 @@ def collect_results(self, write_to_file: bool, snapshot_number: int): self.bs.sinr_ext[bs].tolist(), ) self.results.imt_ul_inr.extend(self.bs.inr[bs].tolist()) + self.results.imt_ul_interf_power.extend( + (self.bs.ext_interference[bs] - 30).tolist(), + ) + self.results.imt_ul_pfd.extend( + self.system_imt_pfd.T[bs].tolist(), + ) self.results.imt_ul_tx_power.extend(self.ue.tx_power[ue].tolist()) imt_ul_tx_power_density = 10 * np.log10( diff --git a/sharc/station_factory.py b/sharc/station_factory.py index 02384d4f2..828435381 100644 --- a/sharc/station_factory.py +++ b/sharc/station_factory.py @@ -10,6 +10,7 @@ import sys import math +from sharc.support.sharc_logger import SimulationLogger from sharc.support.enumerations import StationType from sharc.parameters.parameters import Parameters from sharc.parameters.imt.parameters_imt import ParametersImt @@ -51,6 +52,7 @@ from sharc.antenna.antenna_s1528 import AntennaS1528 from sharc.antenna.antenna_s1855 import AntennaS1855 from sharc.antenna.antenna_s1528 import AntennaS1528, AntennaS1528Leo, AntennaS1528Taylor +from sharc.antenna.antenna_f1245_fs import Antenna_f1245_fs from sharc.antenna.antenna_beamforming_imt import AntennaBeamformingImt from sharc.topology.topology import Topology from sharc.topology.topology_ntn import TopologyNTN @@ -59,6 +61,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 +99,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 +116,9 @@ def generate_imt_base_stations( elev=topology.elevation ) imt_base_stations.is_space_station = True + power_backoff = topology.power_backoff + if param.bs.antenna.pattern == "Antenna System 4": + param.bs.antenna.antenna_system_4.beam_ground_elev_angles = topology.beam_ground_elev_angles else: if topology.determines_local_geometry: imt_base_stations.geom = topology.get_bs_geometry() @@ -132,7 +139,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)], ) @@ -209,6 +218,21 @@ def generate_imt_base_stations( elif param.topology.type == 'HOTSPOT': imt_base_stations.geom.intersite_dist = param.topology.hotspot.intersite_distance + # NOTE: Experimental features for logging and analysis + if param.topology.type == 'MSS_DC': + # número de beams ativos + SimulationLogger.log_to_csv("num_of_active_beams", [np.sum(imt_base_stations.active)]) + # número de satélites, deduplicar posição + all_pos = np.stack((imt_base_stations.geom.x_global, imt_base_stations.geom.y_global, imt_base_stations.geom.z_global), axis=-1) + _, sat_idx, num_of_beams_per_sat = np.unique( + all_pos, + axis=0, + return_index=True, + return_counts=True, + ) + SimulationLogger.log_to_csv("num_of_sat", [len(num_of_beams_per_sat)]) + SimulationLogger.log_to_csv("num_of_beams_per_sat", num_of_beams_per_sat) + return imt_base_stations @staticmethod @@ -869,10 +893,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 +965,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 +1097,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 +1268,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,16 +1330,15 @@ 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": fs_station.antenna = np.array([AntennaOmni(param.antenna_gain)]) elif param.antenna_pattern == "ITU-R F.699": fs_station.antenna = np.array([AntennaF699(param)]) + elif param.antenna_pattern == "ITU-R F.1245_fs": + fs_station.antenna == np.array([Antenna_f1245_fs(param)]) else: sys.stderr.write( "ERROR\nInvalid FS antenna pattern: " + param.antenna_pattern, @@ -1386,6 +1398,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 +1472,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 +1584,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 +1652,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 +1735,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 +1763,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 +1821,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 07ecd3ca6..56fdfd85e 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.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.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.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.pfd_external_aggregated = np.empty(n) # External aggregated PFD in dBm/m^2 + 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/support/sharc_logger.py b/sharc/support/sharc_logger.py index d3f73c9c2..cc7354660 100644 --- a/sharc/support/sharc_logger.py +++ b/sharc/support/sharc_logger.py @@ -1,4 +1,5 @@ # -*- coding: utf-8 -*- +import csv import os import sys import yaml @@ -57,6 +58,41 @@ def get_output_dir(cls) -> Optional[Path]: """Return the global output directory, if set.""" return cls._global_output_dir + @classmethod + def log_to_csv( + cls, + csv_name: str, + vals: list + ): + """ + Log a list of values to a CSV file. + Args: + csv_name (str): The name of the CSV file (without extension). + vals (list): A list of values to be logged to the CSV file. + Returns: + None + Notes: + - Creates a new CSV file if it doesn't exist. + - Appends values to the file if it already exists. + - Automatically writes a "samples" header row for new files. + - Each value is written as a separate row in a single column. + """ + if cls._global_output_dir is None: + # Output directory not yet initialized, skip CSV logging + return + + p = cls._global_output_dir / f"{csv_name}.csv" + + write_header = not p.exists() + + data = [[v] for v in vals] + with open(p, "a", newline="") as file: + writer = csv.writer(file) + if write_header: + writer.writerow(["samples"]) + + writer.writerows(data) + def __init__(self, param_file: str, log_base: str = "simulation_log"): self.param_file: Path = Path(param_file).resolve() self.param_name: str = self.param_file.stem @@ -96,7 +132,7 @@ def end(self): self.output_dir = base_dir / f"simulation_{self.param_name}_{self.timestamp}" self.output_dir.mkdir(parents=True, exist_ok=True) - self.log_path = self.output_dir / f"{self.log_base}_{self.timestamp}.yaml" + self. path = self.output_dir / f"{self.log_base}_{self.timestamp}.yaml" with open(self.log_path, "w") as f: yaml.dump(self.data, f, sort_keys=False, allow_unicode=True) diff --git a/sharc/topology/topology_imt_mss_dc.py b/sharc/topology/topology_imt_mss_dc.py index 337cbe43c..60c8f86cc 100644 --- a/sharc/topology/topology_imt_mss_dc.py +++ b/sharc/topology/topology_imt_mss_dc.py @@ -15,6 +15,7 @@ import geopandas as gpd import functools +from sharc.support.sharc_logger import SimulationLogger from collections import defaultdict from sharc.support.sharc_utils import to_scalar from sharc.topology.topology import Topology @@ -63,6 +64,8 @@ def __init__(self, params: ParametersImtMssDc, self.space_station_x = None self.space_station_y = None self.space_station_z = None + self.beam_ground_elev_angles = None # beam ground elevation angle + self.power_backoff = None # power backoff per beam self.lat = None self.lon = None @@ -278,7 +281,7 @@ def get_coordinates( _, all_azimuth, all_elevation = cartesian_to_polar( pointing_vec_x, pointing_vec_y, pointing_vec_z) - beams_elev, beams_azim, sx, sy = TopologyImtMssDc.get_satellite_pointing( + beams_power_backoff, beams_elev, beams_azim, beams_ground_elev, sx, sy = TopologyImtMssDc.get_satellite_pointing( random_number_gen, coordinate_system, orbit_params, @@ -306,6 +309,18 @@ def get_coordinates( y: list(x) + list(y), beams_azim)) + ground_elev = np.array( + functools.reduce( + lambda x, + y: list(x) + + list(y), + beams_ground_elev)) + space_station_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) @@ -320,6 +335,7 @@ def get_coordinates( 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,)) @@ -341,6 +357,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, @@ -349,6 +366,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) @@ -386,6 +404,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] @@ -411,6 +449,10 @@ def get_satellite_pointing( eligible_sats_msk &= polygon_mask eligible_sats_idx = np.where(eligible_sats_msk)[0] + # NOTE: Experimental features for logging and analysis + # número de satélites candidatos/elegíveis + SimulationLogger.log_to_csv("num_of_candidate_sats", [len(eligible_sats_idx)]) + all_elevations = calc_elevation( grid_lat[:, np.newaxis], all_sat_lat[eligible_sats_msk][np.newaxis, :], @@ -464,23 +506,30 @@ def get_satellite_pointing( # 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( @@ -572,6 +621,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) @@ -581,6 +631,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] @@ -596,7 +650,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( @@ -663,6 +717,12 @@ def calculate_coordinates(self, random_number_gen=np.random.RandomState()): self.y = sat_values["sectors_y"] self.z = sat_values["sectors_z"] + # Beam ground elevation angle + self.beam_ground_elev_angles = sat_values["beams_ground_elev"] + + # 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/e2e/test_integration_imt_victim.py b/tests/e2e/test_integration_imt_victim.py index a7f70e2b1..02700e3fd 100644 --- a/tests/e2e/test_integration_imt_victim.py +++ b/tests/e2e/test_integration_imt_victim.py @@ -141,8 +141,6 @@ def assert_sys_to_imt_dl_results_attr( self.assertEqual(len(res.imt_bs_antenna_gain), n_ue) self.assertEqual(len(res.imt_ue_antenna_gain), n_ue) - self.assertEqual(len(res.imt_dl_pfd_external), n_ue) - self.assertEqual(len(res.imt_dl_pfd_external_aggregated), n_ue) self.assertEqual(len(res.imt_dl_sinr_ext), n_ue) self.assertEqual(len(res.imt_dl_sinr), n_ue) self.assertEqual(len(res.imt_dl_snr), n_ue) @@ -153,6 +151,7 @@ def assert_sys_to_imt_dl_results_attr( # testing attributes that should be per ue towards system self.assertEqual(len(res.imt_system_antenna_gain_adjacent), n_ue * n_sys) self.assertEqual(len(res.imt_system_path_loss), n_ue * n_sys) + self.assertEqual(len(res.imt_dl_pfd), n_ue * n_sys) # NOTE: it may not have co-channel since this test is for adjacent # self.assertEqual(len(res.imt_system_antenna_gain), n_ue * n_sys) @@ -194,6 +193,7 @@ def assert_sys_to_imt_ul_results_attr( self.assertEqual(len(res.imt_ul_inr), n_ue) self.assertEqual(len(res.imt_ul_tput_ext), n_ue) self.assertEqual(len(res.imt_ul_tput), n_ue) + self.assertEqual(len(res.imt_ul_pfd), n_bs * n_sys) # testing attributes that should be per ue towards system self.assertEqual(len(res.imt_system_antenna_gain_adjacent), n_ue * n_sys) diff --git a/tests/parameters/parameters_for_testing.yaml b/tests/parameters/parameters_for_testing.yaml index 86631172b..51f257d1c 100644 --- a/tests/parameters/parameters_for_testing.yaml +++ b/tests/parameters/parameters_for_testing.yaml @@ -283,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 @@ -1024,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", diff --git a/tests/parameters/test_parameters.py b/tests/parameters/test_parameters.py index e5fd967d2..ee614191d 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, @@ -744,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 = [ diff --git a/tests/test_topology_imt_mss_dc.py b/tests/test_topology_imt_mss_dc.py index 69a715ca5..2e9d35f31 100644 --- a/tests/test_topology_imt_mss_dc.py +++ b/tests/test_topology_imt_mss_dc.py @@ -2,7 +2,7 @@ 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 @@ -176,6 +176,9 @@ def test_minimum_service_angle(self): 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(