diff --git a/profile/P619/detailed_profile_p619.py b/profile/P619/detailed_profile_p619.py index 9b2a6ef24..f460bf99b 100644 --- a/profile/P619/detailed_profile_p619.py +++ b/profile/P619/detailed_profile_p619.py @@ -19,10 +19,8 @@ # Constants frequency_MHz = 2680.0 -space_station_alt_m = 35786.0 * 1000 # Geostationary orbit altitude in meters earth_station_alt_m = 1000.0 earth_station_lat_deg = -15.7801 -earth_station_long_diff_deg = 0.0 # Not used in the loss calculation season = "SUMMER" apparent_elevation = np.arange(0, 90) # Elevation angles from 0 to 90 degrees surf_water_vapour_density = 2.5 @@ -30,10 +28,10 @@ # Initialize the propagation model random_number_gen = np.random.RandomState(101) propagation = PropagationP619(random_number_gen=random_number_gen, - space_station_alt_m=space_station_alt_m, earth_station_alt_m=earth_station_alt_m, earth_station_lat_deg=earth_station_lat_deg, - earth_station_long_diff_deg=earth_station_long_diff_deg, + mean_clutter_height='low', + below_rooftop=0.0, season=season) diff --git a/profile/P619/profile_with_and_without_lookuptable_p619.py b/profile/P619/profile_with_and_without_lookuptable_p619.py index 294ecac3c..903e7f8b9 100644 --- a/profile/P619/profile_with_and_without_lookuptable_p619.py +++ b/profile/P619/profile_with_and_without_lookuptable_p619.py @@ -16,7 +16,6 @@ # Constants frequency_MHz = 2680.0 -space_station_alt_m = 35786.0 * 1000 # Geostationary orbit altitude in meters earth_station_alt_m = 1000.0 earth_station_lat_deg = -15.7801 earth_station_long_diff_deg = 0.0 # Not used in the loss calculation @@ -27,10 +26,10 @@ # Initialize the propagation model random_number_gen = np.random.RandomState(101) propagation = PropagationP619(random_number_gen=random_number_gen, - space_station_alt_m=space_station_alt_m, earth_station_alt_m=earth_station_alt_m, earth_station_lat_deg=earth_station_lat_deg, - earth_station_long_diff_deg=earth_station_long_diff_deg, + mean_clutter_height='low', + below_rooftop=0.0, season=season) # Profile with lookup table diff --git a/requirements.txt b/requirements.txt index 039403f36..23bd62a47 100644 --- a/requirements.txt +++ b/requirements.txt @@ -46,14 +46,14 @@ mpmath==1.3.0 multipledispatch==1.0.0 nest-asyncio==1.6.0 numpy==1.26.4 -obspy==1.4.0 +#obspy==1.4.0 overloading==0.5.0 packaging==23.2 -pandas==2.2.1 +pandas==2.3.1 parso==0.8.4 pep8==1.7.1 pexpect==4.9.0 -pillow==10.2.0 +pillow==11.3.0 pip-system-certs==4.0 platformdirs==4.2.0 plotly==5.23.0 @@ -80,9 +80,9 @@ PyYAML==6.0.1 pyzmq==26.0.3 requests==2.31.0 ruff==0.6.9 -scipy==1.12.0 +scipy==1.16.1 setuptools==69.1.1 -shapely==2.0.3 +shapely==2.1.1 six==1.16.0 snowballstemmer==2.2.0 SQLAlchemy==2.0.29 diff --git a/sharc/antenna/antenna_beamforming_imt.py b/sharc/antenna/antenna_beamforming_imt.py index 953ddf923..d041c249e 100644 --- a/sharc/antenna/antenna_beamforming_imt.py +++ b/sharc/antenna/antenna_beamforming_imt.py @@ -14,8 +14,7 @@ from sharc.antenna.antenna_element_imt_const import AntennaElementImtConst from sharc.antenna.antenna_subarray_imt import AntennaSubarrayIMT from sharc.antenna.antenna import Antenna -from sharc.support.named_tuples import AntennaPar -from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt, ParametersAntennaSubarrayImt +from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt class AntennaBeamformingImt(Antenna): @@ -45,11 +44,9 @@ class AntennaBeamformingImt(Antenna): def __init__( self, - par: AntennaPar, + par: ParametersAntennaImt, azimuth: float, - elevation: float, - subarray: ParametersAntennaSubarrayImt = ParametersAntennaSubarrayImt( - is_enabled=False)): + elevation: float): """ Constructs an AntennaBeamformingImt object. Does not receive angles in local coordinate system. @@ -57,7 +54,7 @@ def __init__( Parameters --------- - param (AntennaPar): antenna IMT parameters + param (ParametersAntennaImt): antenna IMT parameters azimuth (float): antenna's physical azimuth inclination elevation (float): antenna's physical elevation inclination referenced in the x axis @@ -78,12 +75,12 @@ def __init__( par.element_pattern} not supported", ) sys.exit(1) - if subarray.is_enabled: + if par.subarray.is_enabled: self.subarray = AntennaSubarrayIMT( element=self.element, - eletrical_downtilt=subarray.eletrical_downtilt, - n_rows=subarray.n_rows, - element_vert_spacing=subarray.element_vert_spacing + eletrical_downtilt=par.subarray.eletrical_downtilt, + n_rows=par.subarray.n_rows, + element_vert_spacing=par.subarray.element_vert_spacing ) self.azimuth = azimuth diff --git a/sharc/antenna/antenna_element_imt_const.py b/sharc/antenna/antenna_element_imt_const.py index 31a88430c..bd9029d5a 100644 --- a/sharc/antenna/antenna_element_imt_const.py +++ b/sharc/antenna/antenna_element_imt_const.py @@ -6,7 +6,7 @@ """ -from sharc.support.named_tuples import AntennaPar +from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt import numpy as np @@ -20,7 +20,7 @@ class AntennaElementImtConst(object): """ - def __init__(self, par: AntennaPar): + def __init__(self, par: ParametersAntennaImt): """ Constructs an AntennaElementImt object. diff --git a/sharc/antenna/antenna_element_imt_f1336.py b/sharc/antenna/antenna_element_imt_f1336.py index 00c00c6a3..443fb046b 100644 --- a/sharc/antenna/antenna_element_imt_f1336.py +++ b/sharc/antenna/antenna_element_imt_f1336.py @@ -8,7 +8,7 @@ import numpy as np import sys -from sharc.support.named_tuples import AntennaPar +from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt class AntennaElementImtF1336(object): @@ -23,7 +23,7 @@ class AntennaElementImtF1336(object): phi_3db (float): horizontal 3dB beamwidth of single element [degrees] """ - def __init__(self, par: AntennaPar): + def __init__(self, par: ParametersAntennaImt): """ Constructs an AntennaElementImt object. diff --git a/sharc/antenna/antenna_element_imt_m2101.py b/sharc/antenna/antenna_element_imt_m2101.py index 3e826e940..f78b7908c 100644 --- a/sharc/antenna/antenna_element_imt_m2101.py +++ b/sharc/antenna/antenna_element_imt_m2101.py @@ -7,7 +7,7 @@ import numpy as np -from sharc.support.named_tuples import AntennaPar +from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt class AntennaElementImtM2101(object): @@ -23,7 +23,7 @@ class AntennaElementImtM2101(object): sla_v (float): element vertical sidelobe attenuation """ - def __init__(self, par: AntennaPar): + def __init__(self, par: ParametersAntennaImt): """ Constructs an AntennaElementImt object. diff --git a/sharc/antenna/antenna_factory.py b/sharc/antenna/antenna_factory.py index e6c0f8b6c..051e1f332 100644 --- a/sharc/antenna/antenna_factory.py +++ b/sharc/antenna/antenna_factory.py @@ -1,9 +1,11 @@ """Antenna factory module for creating antenna instances based on parameters.""" from sharc.parameters.parameters_antenna import ParametersAntenna +from sharc.antenna.antenna import Antenna from sharc.antenna.antenna_mss_adjacent import AntennaMSSAdjacent from sharc.antenna.antenna_omni import AntennaOmni +from sharc.antenna.antenna_mss_hibleo_x_ue import AntennaMssHibleoXUe from sharc.antenna.antenna_f699 import AntennaF699 from sharc.antenna.antenna_s465 import AntennaS465 from sharc.antenna.antenna_rra7_3 import AntennaReg_RR_A7_3 @@ -13,6 +15,8 @@ from sharc.antenna.antenna_s1528 import AntennaS1528, AntennaS1528Leo, AntennaS1528Taylor from sharc.antenna.antenna_beamforming_imt import AntennaBeamformingImt +import numpy as np + class AntennaFactory(): """Factory class for creating antenna instances based on pattern parameters.""" @@ -27,6 +31,8 @@ def create_antenna( match antenna_params.pattern: case "OMNI": return AntennaOmni(antenna_params.gain) + case "HibleoX": + return AntennaMssHibleoXUe(antenna_params.hibleo_x.frequency) case "ITU-R F.699": return AntennaF699(antenna_params.itu_r_f_699) case "ITU-R-S.1528-Taylor": @@ -58,3 +64,35 @@ def create_antenna( raise ValueError( f"Antenna factory does not support pattern { antenna_params.pattern}") + + @staticmethod + def create_n_antennas( + antenna_params: ParametersAntenna, + azimuth: np.ndarray | float, + elevation: np.ndarray | float, + n_stations: int, + ): + """ + Creates many antennas based on passed parameters. + If antenna does not require each object to have different state, + only a single antenna object will be created, and every position + in the array will point to it. + This is much more performant. + """ + antennas = np.empty((n_stations,), dtype=Antenna) + assert n_stations == len(azimuth) + assert n_stations == len(elevation) + + if antenna_params.pattern == "ARRAY": + for i in range(n_stations): + antennas[i] = AntennaFactory.create_antenna( + antenna_params, azimuth[i], elevation[i], + ) + else: + # some antennas don't need azimuth and elevation at all + # this makes it much faster + antennas[:] = AntennaFactory.create_antenna( + antenna_params, None, None, + ) + + return antennas diff --git a/sharc/antenna/antenna_mss_adjacent.py b/sharc/antenna/antenna_mss_adjacent.py index b624e2ebe..840a4f928 100644 --- a/sharc/antenna/antenna_mss_adjacent.py +++ b/sharc/antenna/antenna_mss_adjacent.py @@ -43,9 +43,10 @@ def calculate_gain(self, *args, **kwargs) -> np.array: np.array Calculated antenna gain values. """ - theta = np.absolute(kwargs["off_axis_angle_vec"]) + theta_rad = np.deg2rad(np.absolute(kwargs["off_axis_angle_vec"])) + theta_rad = np.minimum(theta_rad, np.pi / 2 - 1e-4) return 20 * np.log10(self.frequency / 2e3) + 10 * \ - np.log10(np.cos(np.deg2rad(theta))) + np.log10(np.cos(theta_rad)) if __name__ == '__main__': diff --git a/sharc/antenna/antenna_mss_hibleo_x_ue.py b/sharc/antenna/antenna_mss_hibleo_x_ue.py new file mode 100644 index 000000000..10efaae32 --- /dev/null +++ b/sharc/antenna/antenna_mss_hibleo_x_ue.py @@ -0,0 +1,128 @@ + +# -*- coding: utf-8 -*- +"""Antenna model for HibleoX MSS UE in the 2.5GHz frequency band.""" +from sharc.antenna.antenna import Antenna + +import numpy as np + +# The filing does not provide an antenna model, but a set of points. +# We interpolate these points to create a continuous function. +# The points are taken from the filing. +data_str = """ +x, y +0, 2.7 +10, 2.6566311598597308 +20, 2.431925804835606 +30, 2.150602473172147 +40, 1.6546431302238513 +50, 0.9929370896722656 +60, 0.0904942389326866 +70, -0.9031468347087825 +80, -2.386785139873968 +90, -3.679459224299312 +""" + +HALF_BEAM_WIDTH_DEG = 64.0 + + +class AntennaMssHibleoXUe(Antenna): + """ + Implements the MSS UE antenna pattern provided for the HIBLEO-X system in it's filings. + + This is provided as a reference antenna pattern in document Document 4C/166-E - 30 September 2024. + This is valid only for frequencies around 2.45 GHz. + """ + + def __init__(self, frequency_MHz: float): + """ + Initialize the AntennaMSSAdjacent class. + + Parameters + ---------- + frequency_MHz : float + Frequency in MHz for the antenna model. + """ + super().__init__() + if frequency_MHz < 2483 or frequency_MHz > 2500: + raise ValueError("This antenna model is only valid for frequencies around 2.45 GHz") + self.frequency = frequency_MHz + + # Parse the data string to extract x and y values + lines = data_str.strip().split('\n')[1:] + self.pattern_theta_deg = np.array([float(line.split(',')[0]) for line in lines]) + self.patten_gain_dbi = np.array([float(line.split(',')[1]) for line in lines]) + + def calculate_gain(self, *args, **kwargs) -> np.array: + """ + Calculate the antenna gain for the given off-axis angles. + + Parameters + ---------- + *args : tuple + Positional arguments (not used). + **kwargs : dict + Keyword arguments, expects 'off_axis_angle_vec' in degrees as input. + + Returns + ------- + np.array + Calculated antenna gain values. + """ + theta_deg = kwargs["off_axis_angle_vec"] + return -2.5226404 * (np.deg2rad(theta_deg) ** 2) + 2.7 + + +if __name__ == '__main__': + import plotly.graph_objects as go + from scipy.optimize import curve_fit + + # Compare the fitted parabolic function to the provided points + + mss_ue_antenna = AntennaMssHibleoXUe(2483) + + # Define the custom parabolic function without the linear term + def poly_no_linear(x, a): + """ + Polynomial function of the form ax^2 + c (no linear term). + """ + return a * x**2 + np.max(mss_ue_antenna.patten_gain_dbi) + + # Fit the data using curve_fit + # initial_guess = [1, 1] # Optional: provide an initial guess for the parameters + params, covariance = curve_fit( + poly_no_linear, + np.deg2rad(mss_ue_antenna.pattern_theta_deg), + mss_ue_antenna.patten_gain_dbi + ) + + print("Fitting parabolic function (no linear term) to the provided points:") + print(f"Fitted parameters:\na={params}") + print(f"Covariance matrix:\n{covariance}") + + off_axis_vec_deg = np.linspace(-90, 90, 100) + gains = mss_ue_antenna.calculate_gain(off_axis_angle_vec=off_axis_vec_deg) + + gains_interp = np.interp( + np.abs(off_axis_vec_deg), + mss_ue_antenna.pattern_theta_deg, + mss_ue_antenna.patten_gain_dbi + ) + + fig = go.Figure(data=go.Scatter(x=off_axis_vec_deg, y=gains, mode='markers+lines', name='Parabolic Fit pattern')) + fig.update_layout( + title='Hibleo-X UE pattern from filing', + xaxis_title='off-axis angle (degrees)', + yaxis_title='Gain (dBi)', + template="plotly_white" + ) + fig.update_yaxes( + linewidth=1, + linecolor='black', + mirror=True, + ticks='inside', + showline=True, + gridcolor="#DCDCDC", + gridwidth=1.5, + ) + fig.add_trace(go.Scatter(x=off_axis_vec_deg, y=gains_interp, mode='markers', name='Interpolated')) + fig.show() diff --git a/sharc/antenna/antenna_subarray_imt.py b/sharc/antenna/antenna_subarray_imt.py index eaaffb623..0b8afc8ae 100644 --- a/sharc/antenna/antenna_subarray_imt.py +++ b/sharc/antenna/antenna_subarray_imt.py @@ -206,7 +206,7 @@ def calculate_gain( # Plot BS TX radiation patterns par = bs_param.get_antenna_parameters() - bs_array = AntennaBeamformingImt(par, 0, 0, bs_param.subarray) + bs_array = AntennaBeamformingImt(par, 0, -6) f = plot.plot_element_pattern(bs_array, "BS", "ELEMENT") # f.savefig(figs_dir + "BS_element.pdf", bbox_inches='tight') f = plot.plot_element_pattern(bs_array, "BS", "SUBARRAY") @@ -215,7 +215,7 @@ def calculate_gain( # Plot UE TX radiation patterns par = bs2_param.get_antenna_parameters() - bs2_array = AntennaBeamformingImt(par, 0, 0, bs2_param.subarray) + bs2_array = AntennaBeamformingImt(par, 0, 0) # plot.plot_element_pattern(bs2_array, "BS 2", "ELEMENT") # plot.plot_element_pattern(bs2_array, "BS 2", "ARRAY") diff --git a/sharc/antenna/beamforming_normalization/beamforming_normalizer.py b/sharc/antenna/beamforming_normalization/beamforming_normalizer.py index c99498e59..02aaaf57b 100644 --- a/sharc/antenna/beamforming_normalization/beamforming_normalizer.py +++ b/sharc/antenna/beamforming_normalization/beamforming_normalizer.py @@ -10,7 +10,7 @@ from sys import stdout from sharc.antenna.antenna_beamforming_imt import AntennaBeamformingImt -from sharc.support.named_tuples import AntennaPar +from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt class BeamformingNormalizer(object): @@ -74,7 +74,7 @@ def __init__(self, res_deg: float, tol: float): def generate_correction_matrix( self, - par: AntennaPar, + par: ParametersAntennaImt, file_name: str, testing=False, ): @@ -82,7 +82,7 @@ def generate_correction_matrix( Generates the correction factor matrix and saves it in a file Parameters: - par (AntennaPar): set of antenna parameters to which calculate the + par (ParametersAntennaImt): set of antenna parameters to which calculate the correction factor file_name (str): name of file to which save the correction matrix """ @@ -200,7 +200,7 @@ def _save_files(self, cf_co, err_co, cf_adj, err_adj, par, file_name): antenna element error_adj_channel (tuple): lower and upper bounds [dB] of single antenna element correction factor - parameters (AntennaPar): antenna parameters used in the + parameters (ParametersAntennaImt): antenna parameters used in the normalization Parameters: @@ -239,35 +239,22 @@ def _save_files(self, cf_co, err_co, cf_adj, err_adj, par, file_name): norm = BeamformingNormalizer(resolution, tolerance) # Antenna parameters - adjacent_antenna_model = "SINGLE_ELEMENT" - normalization = False - norm_file = None - element_pattern = "M2101" - element_max_g = 5 - element_phi_deg_3db = 65 - element_theta_deg_3db = 65 - element_am = 30 - element_sla_v = 30 - n_rows = 8 - n_columns = 8 - horiz_spacing = 0.5 - vert_spacing = 0.5 - minimum_array_gain = -200 - down_tilt = 0 - par = AntennaPar( - normalization, - norm_file, - element_pattern, - element_max_g, - element_phi_deg_3db, - element_theta_deg_3db, - element_am, - element_sla_v, - n_rows, - n_columns, - horiz_spacing, - vert_spacing, - down_tilt, + par = ParametersAntennaImt( + adjacent_antenna_model="SINGLE_ELEMENT", + normalization=False, + norm_file=None, + element_pattern="M2101", + element_max_g=5, + element_phi_3db=65, + element_theta_3db=65, + element_am=30, + element_sla_v=30, + n_rows=8, + n_columns=8, + element_horiz_spacing=0.5, + element_vert_spacing=0.5, + minimum_array_gain=-200, + downtilt=0, ) # Set range of values & calculate correction factor diff --git a/sharc/antenna/beamforming_normalization/normalize_script.py b/sharc/antenna/beamforming_normalization/normalize_script.py index 71084570c..58fd6872f 100644 --- a/sharc/antenna/beamforming_normalization/normalize_script.py +++ b/sharc/antenna/beamforming_normalization/normalize_script.py @@ -25,7 +25,7 @@ norm (BeamformingNormalizer): object that calculates the normalization. param_list (list): list of antenna parameters to which calculate the correction factors. New parameters are added as: - AntennaPar(adjacent_antenna_model, + ParametersAntennaImt(adjacent_antenna_model, normalization, norm_data, element_pattern, @@ -43,7 +43,7 @@ normalization parameter must be set to False, otherwise script will try to normalize an already normalized antenna. file_names (list): list of file names to which save the normalization data. - Files are paired with AntennaPar objects in param_list, so that the + Files are paired with ParametersAntennaImt objects in param_list, so that the normalization data of the first element of param_list is saved in a file with the name specified in the first element of file_names and so on. @@ -64,11 +64,11 @@ element error_adj_channel (tuple): lower and upper bounds [dB] of single antenna element correction factor - parameters (AntennaPar): antenna parameters used in the normalization + parameters (ParametersAntennaImt): antenna parameters used in the normalization """ from sharc.antenna.beamforming_normalization.beamforming_normalizer import BeamformingNormalizer -from sharc.support.named_tuples import AntennaPar +from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt import sys import os sys.path.append(os.path.join(os.path.dirname(__file__), "..", "..", "..")) @@ -78,42 +78,28 @@ ########################################################################### # List of antenna parameters to which calculate the normalization factors. - adjacent_antenna_model = "" # not needed here - normalization = False # not needed here - normalization_data = None # not needed here - element_pattern = "M2101" - element_max_g = 5 - element_phi_3db = 65 - element_theta_3db = 65 - element_am = 30 - element_sla_v = 30 - n_rows = 8 - n_columns = 8 - element_horiz_spacing = 0.5 - element_vert_spacing = 0.5 - multiplication_factor = 12 - minimum_array_gain = -200 - downtilt = 0 - file_names = ["bs_norm_8x8_050.npz"] param_list = [ - AntennaPar( - adjacent_antenna_model, - normalization, - normalization_data, - element_pattern, - element_max_g, - element_phi_3db, - element_theta_3db, - element_am, - element_sla_v, - n_rows, - n_columns, - element_horiz_spacing, - element_vert_spacing, - multiplication_factor, - minimum_array_gain, - downtilt, + ParametersAntennaImt( + # not needed here + adjacent_antenna_model="", + # not needed here + normalization=False, + # not needed here + normalization_file=None, + element_pattern="M2101", + element_max_g=5, + element_phi_3db=65, + element_theta_3db=65, + element_am=30, + element_sla_v=30, + n_rows=8, + n_columns=8, + element_horiz_spacing=0.5, + element_vert_spacing=0.5, + multiplication_factor=12, + minimum_array_gain=-200, + downtilt=0, ), ] ########################################################################### diff --git a/sharc/campaigns/mss_d2d_to_imt/scripts/orbit_model_anaylisis.py b/sharc/campaigns/mss_d2d_to_imt/scripts/orbit_model_anaylisis.py index f264adc1c..38c6e504c 100644 --- a/sharc/campaigns/mss_d2d_to_imt/scripts/orbit_model_anaylisis.py +++ b/sharc/campaigns/mss_d2d_to_imt/scripts/orbit_model_anaylisis.py @@ -48,7 +48,10 @@ delta=orbit_params.inclination_deg, hp=orbit_params.perigee_alt_km, ha=orbit_params.apogee_alt_km, - Mo=orbit_params.initial_mean_anomaly + Mo=orbit_params.initial_mean_anomaly, + model_time_as_random_variable=False, + t_min=0.0, + t_max=None ) # Show visible satellites from ground-station diff --git a/sharc/input/parameters.yaml b/sharc/input/parameters.yaml index f9c384041..8edf9877f 100644 --- a/sharc/input/parameters.yaml +++ b/sharc/input/parameters.yaml @@ -292,6 +292,9 @@ imt: sats_per_plane: 32 long_asc_deg: 18.0 phasing_deg: 3.9 + model_time_as_random_variable: False + t_min: 0.0 + t_max: null - n_planes: 12 inclination_deg: 26.0 perigee_alt_km: 580.0 @@ -299,6 +302,9 @@ imt: sats_per_plane: 20 long_asc_deg: 30.0 phasing_deg: 2.0 + model_time_as_random_variable: False + t_min: 0.0 + t_max: null - n_planes: 26 inclination_deg: 97.77 perigee_alt_km: 595 @@ -306,6 +312,9 @@ imt: sats_per_plane: 30 long_asc_deg: 14.0 phasing_deg: 7.8 + model_time_as_random_variable: False + t_min: 0.0 + t_max: null ########################################################################### # Defines the antenna model to be used in compatibility studies between # IMT and other services in adjacent band @@ -606,6 +615,11 @@ imt: # difference between longitudes of IMT-NTN station and FSS-ES system [deg] # (positive if space-station is to the East of earth-station) earth_station_long_diff_deg: 0 + # Average height of clutter. According to 2018 it can be "Low", "Mid" or "High" + mean_clutter_height: "High" + # Percentage of stations below rooftop (To apply clutter) + below_rooftop: 100 + single_earth_station: ########################################################################### # Sensor center frequency [MHz] @@ -751,6 +765,10 @@ single_earth_station: # difference between longitudes of IMT-NTN station and FSS-ES system [deg] # (positive if space-station is to the East of earth-station) earth_station_long_diff_deg: 10 + # Average height of clutter. According to 2018 it can be "Low", "Mid" or "High" + mean_clutter_height: "High" + # Percentage of stations below rooftop (To apply clutter) + below_rooftop: 100 ########################################################################### # season - season of the year. @@ -799,7 +817,8 @@ single_earth_station: ########################################################################### # determine whether clutter loss following ITU-R P.2108 is added (TRUE/FALSE) clutter_loss: TRUE - + # Determine if clutter is applied to "one_end" or "both)ends" + clutter_type: "one_end" # HDFSS propagation parameters param_hdfss: ########################################################################### @@ -1097,7 +1116,10 @@ mss_d2d: ########################################################################### # year season: SUMMER of WINTER season: SUMMER - + # Average height of clutter. According to 2018 it can be "Low", "Mid" or "High" + mean_clutter_height: "High" + # Percentage of stations below rooftop (To apply clutter) + below_rooftop: 100 # Orbit parameters orbit: # Number of planes diff --git a/sharc/mask/spectral_mask_mss.py b/sharc/mask/spectral_mask_mss.py index c0891278d..deebcb3d8 100644 --- a/sharc/mask/spectral_mask_mss.py +++ b/sharc/mask/spectral_mask_mss.py @@ -79,12 +79,12 @@ def set_mask(self, p_tx): Set the spectral mask (mask_dbm attribute) based on station type, operating frequency and transmit power. Parameters: - p_tx (float): station transmit power. + p_tx (float): station transmit power in dBm """ # dBm/MHz # this should work for the document's dBsd definition # when we have a uniform PSD in assigned band - self.p_tx = p_tx - 10 * np.log10(self.band_mhz) + 30 + self.p_tx = p_tx - 10 * np.log10(self.band_mhz) # attenuation mask mask_dbsd = 40 * np.log10( diff --git a/sharc/parameters/imt/parameters_antenna_imt.py b/sharc/parameters/imt/parameters_antenna_imt.py index c846053ef..a2eed50d3 100644 --- a/sharc/parameters/imt/parameters_antenna_imt.py +++ b/sharc/parameters/imt/parameters_antenna_imt.py @@ -5,7 +5,6 @@ @author: Calil """ -from sharc.support.named_tuples import AntennaPar from numpy import load import typing @@ -58,8 +57,8 @@ class ParametersAntennaImt(ParametersBase): # PS: it isn't implemented for UEs # and current implementation doesn't make sense for UEs horizontal_beamsteering_range: tuple[float | - int, float | int] = (-180., 180.) - vertical_beamsteering_range: tuple[float | int, float | int] = (0., 180.) + int, float | int] = (-180., 179.9999) + vertical_beamsteering_range: tuple[float | int, float | int] = (0., 179.9999) # Mechanical downtilt [degrees]. # PS: downtilt doesn't make sense on UE's @@ -96,6 +95,9 @@ class ParametersAntennaImt(ParametersBase): subarray: ParametersAntennaSubarrayImt = field( default_factory=ParametersAntennaSubarrayImt) + def __post_init__(self): + self.normalization_data = None + def load_subparameters(self, ctx: str, params: dict, quiet=True): """ Load parameters when this class is used as a subparameter. @@ -172,11 +174,11 @@ def validate(self, ctx: str): f"Invalid {ctx}.horizontal_beamsteering_range={self.horizontal_beamsteering_range}\n" "The second value must be bigger than the first" ) - if not all(map(lambda x: x >= -180. and x <= 180., + if not all(map(lambda x: x >= -180. and x < 180., self.horizontal_beamsteering_range)): raise ValueError( f"Invalid {ctx}.horizontal_beamsteering_range={self.horizontal_beamsteering_range}\n" - "Horizontal beamsteering limit angles must be in the range [-180, 180]" + "Horizontal beamsteering limit angles must be in the range [-180, 180)" ) if isinstance(self.vertical_beamsteering_range, list): @@ -197,21 +199,16 @@ def validate(self, ctx: str): f"Invalid {ctx}.vertical_beamsteering_range={self.vertical_beamsteering_range}\n" "The second value must be bigger than the first" ) - if not all(map(lambda x: x >= 0. and x <= 180., + if not all(map(lambda x: x >= 0. and x < 180., self.vertical_beamsteering_range)): raise ValueError( f"Invalid {ctx}.vertical_beamsteering_range={self.vertical_beamsteering_range}\n" - "vertical beamsteering limit angles must be in the range [0, 180]" + "vertical beamsteering limit angles must be in the range [0, 180)" ) - def get_antenna_parameters(self) -> AntennaPar: + def get_normalization_data_if_needed(self): """ - Get the antenna parameters as an AntennaPar object. - - Returns - ------- - AntennaPar - The antenna parameters object constructed from the current configuration. + This loads normalization data if normalization should be applied """ if self.normalization: # Load data, save it in dict and close it @@ -221,23 +218,16 @@ def get_antenna_parameters(self) -> AntennaPar: data.close() else: self.normalization_data = None - tpl = AntennaPar( - self.adjacent_antenna_model, - self.normalization, - self.normalization_data, - self.element_pattern, - self.element_max_g, - self.element_phi_3db, - self.element_theta_3db, - self.element_am, - self.element_sla_v, - self.n_rows, - self.n_columns, - self.element_horiz_spacing, - self.element_vert_spacing, - self.multiplication_factor, - self.minimum_array_gain, - self.downtilt, - ) - - return tpl + + def get_antenna_parameters(self) -> "ParametersAntennaImt": + """ + Get the antenna parameters loadind normalization values if needed. + + Returns + ------- + ParametersAntennaImt + The antenna parameters object constructed from the current configuration. + """ + self.get_normalization_data_if_needed() + + return self diff --git a/sharc/parameters/imt/parameters_imt.py b/sharc/parameters/imt/parameters_imt.py index 5fd077c78..1d007b540 100644 --- a/sharc/parameters/imt/parameters_imt.py +++ b/sharc/parameters/imt/parameters_imt.py @@ -110,11 +110,11 @@ class ParametersUE(ParametersBase): def validate(self, ctx: str): """Validate the UE antenna beamsteering range parameters.""" - if self.antenna.array.horizontal_beamsteering_range != (-180., 180.)\ - or self.antenna.array.vertical_beamsteering_range != (0., 180.): + if self.antenna.array.horizontal_beamsteering_range != (-180., 179.9999)\ + or self.antenna.array.vertical_beamsteering_range != (0., 179.9999): raise NotImplementedError( "UE antenna beamsteering limit has not been implemented. Default values of\n" - "horizontal = (-180., 180.), vertical = (0., 180.) should not be changed") + "horizontal = (-180., 179.9999), vertical = (0., 179.9999) should not be changed") ue: ParametersUE = field(default_factory=ParametersUE) @@ -144,14 +144,10 @@ class ParamatersDL(ParametersBase): # the BS space station and the UEs on Earth's surface. # For now, the NTN footprint is centered over the BS nadir point, therefore # the paramters imt_lag_deg and imt_long_diff_deg SHALL be zero. - # space_station_alt_m - altitude of IMT space station (meters) # earth_station_alt_m - altitude of IMT earth stations (UEs) (in meters) # earth_station_lat_deg - latitude of IMT earth stations (UEs) (in degrees) - # earth_station_long_diff_deg - difference between longitudes of IMT space and earth stations - # (positive if space-station is to the East of earth-station) # season - season of the year. param_p619 = ParametersP619() - season: str = "SUMMER" # TODO: create parameters for where this is needed los_adjustment_factor: float = 18.0 @@ -201,11 +197,6 @@ def load_parameters_from_file(self, config_file: str): f"ParametersImt: Invalid channel model { self.channel_model} for topology NTN", ) - if self.season not in ["SUMMER", "WINTER"]: - raise ValueError(f"ParamtersImt: \ - Invalid value for parameter season - {self.season}. \ - Possible values are \"SUMMER\", \"WINTER\".") - if self.topology.type == "NTN": self.is_space_to_earth = True # self.param_p619.load_from_paramters(self) diff --git a/sharc/parameters/imt/parameters_imt_mss_dc.py b/sharc/parameters/imt/parameters_imt_mss_dc.py index cf69f86eb..3f31ed8da 100644 --- a/sharc/parameters/imt/parameters_imt_mss_dc.py +++ b/sharc/parameters/imt/parameters_imt_mss_dc.py @@ -136,6 +136,8 @@ class ParametersServiceGrid(ParametersBase): country_names: list[str] = field(default_factory=lambda: list([""])) + transform_grid_randomly: bool = False + # margin from inside of border [km] # if positive, makes border smaller by x km # if negative, makes border bigger by x km @@ -197,7 +199,7 @@ def validate(self, ctx: str): raise ValueError( f"{ctx}.eligible_sats_margin_from_border needs to be a number") - self.reset_grid(ctx) + self._load_geom_from_file_if_needed(ctx) super().validate(ctx) @@ -218,13 +220,30 @@ def load_from_active_sat_conditions( if self.eligible_sats_margin_from_border is None: self.eligible_sats_margin_from_border = sat_is_active_if.lat_long_inside_country.margin_from_border - def reset_grid(self, ctx: str, force_update=False): + def reset_grid( + self, + ctx: str, + rng: np.random.RandomState, + force_update=False, + ): """ After creating grid, there are some features that can only be implemented - with knowledge of other parts of the simulator. This method's purpose is - to run only once at the start of the simulation + with knowledge of other parts of the simulator. """ - if self.lon_lat_grid is not None and not force_update: + self._load_geom_from_file_if_needed(ctx, force_update) + + self.lon_lat_grid = generate_grid_in_multipolygon( + self.grid_borders_polygon, + self.beam_radius, + self.transform_grid_randomly, + rng + ) + + self.ecef_grid = lla2ecef( + self.lon_lat_grid[1], self.lon_lat_grid[0], 0) + + def _load_geom_from_file_if_needed(self, ctx: str, force_update=False): + if self.eligibility_polygon is not None and not force_update: return filtered_gdf = load_gdf( self.country_shapes_filename, @@ -239,17 +258,13 @@ def reset_grid(self, ctx: str, force_update=False): shrinked = shrink_countries_by_km( filtered_gdf.geometry.values, self.grid_margin_from_border ) - polygon = shp.ops.unary_union(shrinked) - assert polygon.is_valid, shp.validation.explain_validity(polygon) - assert not polygon.is_empty, "Can't have a empty polygon as filter" + self.grid_borders_polygon = shp.ops.unary_union(shrinked) - self.lon_lat_grid = generate_grid_in_multipolygon( - polygon, - self.beam_radius - ) + assert self.grid_borders_polygon.is_valid, \ + shp.validation.explain_validity(self.grid_borders_polygon) - self.ecef_grid = lla2ecef( - self.lon_lat_grid[1], self.lon_lat_grid[0], 0) + assert not self.grid_borders_polygon.is_empty, \ + "Can't have a empty grid_borders_polygon as filter" self.eligibility_polygon = shp.ops.unary_union(shrink_countries_by_km( filtered_gdf.geometry.values, self.eligible_sats_margin_from_border diff --git a/sharc/parameters/imt/parameters_single_bs.py b/sharc/parameters/imt/parameters_single_bs.py index 4feed03cd..00c970d88 100644 --- a/sharc/parameters/imt/parameters_single_bs.py +++ b/sharc/parameters/imt/parameters_single_bs.py @@ -1,4 +1,5 @@ from dataclasses import dataclass +from typing import List from sharc.parameters.parameters_base import ParametersBase @@ -11,6 +12,7 @@ class ParametersSingleBS(ParametersBase): intersite_distance: int = None cell_radius: int = None num_clusters: int = 1 + azimuth: List[float] | str | None = None def load_subparameters(self, ctx: str, params: dict, quiet=True): """ @@ -61,6 +63,19 @@ def validate(self, ctx): if self.num_clusters not in [1, 2]: raise ValueError(f"{ctx}.num_clusters should either be 1 or 2") + if self.azimuth is None: + self.azimuth = [0.0, 180.0][:self.num_clusters] + elif isinstance(self.azimuth, list): + if len(self.azimuth) != self.num_clusters: + raise ValueError(f"{ctx}.azimuth length must be equal to num_clusters") + if not all(isinstance(a, float) for a in self.azimuth): + raise TypeError(f"{ctx}.azimuth must be a list of floats or a \"random\"") + elif isinstance(self.azimuth, str): + if self.azimuth.lower() != "random": + raise ValueError(f"{ctx}.azimuth must be a list of floats or a \"random\"") + else: + raise TypeError(f"{ctx}.azimuth must be a list of floats or a \"random\"") + if __name__ == "__main__": # Run validation for input parameters diff --git a/sharc/parameters/parameters_antenna.py b/sharc/parameters/parameters_antenna.py index 34e16cda0..f004b3874 100644 --- a/sharc/parameters/parameters_antenna.py +++ b/sharc/parameters/parameters_antenna.py @@ -17,6 +17,7 @@ class ParametersAntenna(ParametersBase): # available antenna radiation patterns __SUPPORTED_ANTENNA_PATTERNS = [ "OMNI", + "HibleoX", "ITU-R F.699", "ITU-R S.465", "ITU-R S.580", @@ -31,6 +32,7 @@ class ParametersAntenna(ParametersBase): # chosen antenna radiation pattern pattern: typing.Literal["OMNI", + "HibleoX", "ITU-R F.699", "ITU-R S.465", "ITU-R S.580", @@ -50,6 +52,10 @@ class ParametersAntenna(ParametersBase): default_factory=ParametersAntennaWithFreq, ) + hibleo_x: ParametersAntennaWithFreq = field( + default_factory=ParametersAntennaWithFreq, + ) + itu_r_f_699: ParametersAntennaWithDiameter = field( default_factory=ParametersAntennaWithDiameter, ) @@ -151,6 +157,8 @@ def validate(self, ctx): match self.pattern: case "OMNI": pass + case "HibleoX": + self.hibleo_x.validate(f"{ctx}.hibleo_x") case "ITU-R F.699": self.itu_r_f_699.validate(f"{ctx}.itu_r_f_699") case "ITU-R S.465": diff --git a/sharc/parameters/parameters_fss_es.py b/sharc/parameters/parameters_fss_es.py index 626bbf759..50aa5fbb8 100644 --- a/sharc/parameters/parameters_fss_es.py +++ b/sharc/parameters/parameters_fss_es.py @@ -97,19 +97,15 @@ class ParametersFssEs(ParametersBase): # Determine whether clutter loss following ITU-R P.2108 is added # (TRUE/FALSE) clutter_loss: bool = True + clutter_type: str = "one-end" # Parameters for the P.619 propagation model used for sharing studies between IMT-NTN and FSS-ES - # space_station_alt_m - altiteude of the IMT-MSS station # earth_station_alt_m - altitude of FSS-ES system (in meters) # earth_station_lat_deg - latitude of FSS-ES system (in degrees) - # earth_station_long_diff_deg - difference between longitudes of IMT-NTN station and FSS-ES system - # (positive if space-station is to the East of earth-station) # season - season of the year. param_p619 = ParametersP619() - space_station_alt_m: float = 35780000.0 earth_station_alt_m: float = 0.0 earth_station_lat_deg: float = 0.0 - earth_station_long_diff_deg: float = 0.0 season: str = "SUMMER" # HDFSS propagation parameters diff --git a/sharc/parameters/parameters_fss_ss.py b/sharc/parameters/parameters_fss_ss.py index 8e5acda0f..434a5205c 100644 --- a/sharc/parameters/parameters_fss_ss.py +++ b/sharc/parameters/parameters_fss_ss.py @@ -29,6 +29,16 @@ class ParametersFssSs(ParametersBase): elevation: float = 270.0 # Azimuth angle [deg] azimuth: float = 0.0 + # Parameters for geometry of the space station + # Space station altitude [m] + space_station_alt_m: float = 35780000.0 + # Earth station altitude [m] + earth_station_alt_m: float = 0.0 + # Earth statin latitude [deg] + earth_station_lat_deg: float = 0.0 + # Difference between longituted of the space station and the earth station [deg] + # (positive if space-station is to the East of earth-station) + earth_station_long_diff_deg: float = 0.0 # System receive noise temperature [K] noise_temperature: float = 950.0 # Adjacent channel reception type. @@ -52,16 +62,8 @@ class ParametersFssSs(ParametersBase): ############################ # Parameters for the P.619 propagation model - # earth_station_alt_m - altitude of IMT system (in meters) - # earth_station_lat_deg - latitude of IMT system (in degrees) - # earth_station_long_diff_deg - difference between longitudes of IMT and satellite system - # (positive if space-station is to the East of earth-station) - # season - season of the year. + param_p619 = ParametersP619() - space_station_alt_m: float = 35780000.0 - earth_station_alt_m: float = 0.0 - earth_station_lat_deg: float = 0.0 - earth_station_long_diff_deg: float = 0.0 season: str = "SUMMER" # Channel parameters # channel model, possible values are "FSPL" (free-space path loss), @@ -104,7 +106,11 @@ def load_parameters_from_file(self, config_file: str): raise ValueError(f"ParametersFssSs: \ Invalid value for paramter channel_model = {self.channel_model}. \ Possible values are \"FSPL\", \"SatelliteSimple\", \"P619\".") - self.param_p619.load_from_paramters(self) + + allowed = {"low", "mid", "high"} + if str(self.param_p619.mean_clutter_height).lower() not in allowed: + raise ValueError("Invalid type of mean_clutter_height. mean_clutter_height must be 'Low', 'Mid', or 'High'") + self.antenna_s1528.set_external_parameters( frequency=self.frequency, bandwidth=self.bandwidth, diff --git a/sharc/parameters/parameters_mss_d2d.py b/sharc/parameters/parameters_mss_d2d.py index 7b5c7608a..172964dbb 100644 --- a/sharc/parameters/parameters_mss_d2d.py +++ b/sharc/parameters/parameters_mss_d2d.py @@ -94,10 +94,6 @@ class ParametersMssD2d(ParametersBase): # paramters for channel model param_p619: ParametersP619 = field(default_factory=ParametersP619) - earth_station_alt_m: float = 0.0 - earth_station_lat_deg: float = 0.0 - earth_station_long_diff_deg: float = 0.0 - season: str = "SUMMER" # Channel parameters # channel model, possible values are "FSPL" (free-space path loss), # "SatelliteSimple" (FSPL + 4 + clutter loss) @@ -200,21 +196,6 @@ def propagate_parameters(self): self.sat_is_active_if, ) - if self.channel_model == "P619": - # mean station altitude in meters - m_alt = 0 - for orbit in self.orbits: - m_alt += orbit.perigee_alt_km * 1e3 - m_alt /= len(self.orbits) - - self.param_p619.set_external_parameters( - space_station_alt_m=m_alt, - earth_station_alt_m=self.earth_station_alt_m, - earth_station_lat_deg=self.earth_station_lat_deg, - earth_station_long_diff_deg=self.earth_station_long_diff_deg, - season=self.season - ) - if __name__ == "__main__": # Run validation for input parameters diff --git a/sharc/parameters/parameters_mss_ss.py b/sharc/parameters/parameters_mss_ss.py index 3419ffb69..1a2b7d668 100644 --- a/sharc/parameters/parameters_mss_ss.py +++ b/sharc/parameters/parameters_mss_ss.py @@ -81,10 +81,7 @@ class ParametersMssSs(ParametersBase): # paramters for channel model param_p619: ParametersP619 = field(default_factory=ParametersP619) - space_station_alt_m: float = 35780000.0 - earth_station_alt_m: float = 0.0 - earth_station_lat_deg: float = 0.0 - earth_station_long_diff_deg: float = 0.0 + season: str = "SUMMER" # Channel parameters # channel model, possible values are "FSPL" (free-space path loss), @@ -147,12 +144,3 @@ def load_parameters_from_file(self, config_file: str): raise ValueError( f"Invalid channel model name { self.channel_model}") - - if self.channel_model == "P619": - self.param_p619.set_external_parameters( - space_station_alt_m=self.altitude, - earth_station_alt_m=self.earth_station_alt_m, - earth_station_lat_deg=self.earth_station_lat_deg, - earth_station_long_diff_deg=self.earth_station_long_diff_deg, - season=self.season - ) diff --git a/sharc/parameters/parameters_ngso_constellation.py b/sharc/parameters/parameters_ngso_constellation.py index b0a08852b..7b2a263c7 100644 --- a/sharc/parameters/parameters_ngso_constellation.py +++ b/sharc/parameters/parameters_ngso_constellation.py @@ -56,13 +56,6 @@ class ParametersNgsoConstellation(ParametersBase): # Channel model configuration # Parameters for the P.619 channel model param_p619: ParametersP619 = field(default_factory=ParametersP619) - # Altitude of the space station in meters - space_station_alt_m: float = 35780000.0 - earth_station_alt_m: float = 0.0 # Altitude of the earth station in meters - earth_station_lat_deg: float = 0.0 # Latitude of the earth station in degrees - # Longitude difference for the earth station in degrees - earth_station_long_diff_deg: float = 0.0 - season: str = "SUMMER" # Season for atmospheric conditions # Channel model type (e.g., FSPL, SatelliteSimple, P619) channel_model: str = "P619" @@ -131,15 +124,6 @@ def load_parameters_from_file(self, config_file: str): f"ParametersNgsoMss: Invalid channel model name { self.channel_model}") - if self.channel_model == "P619": - self.param_p619.set_external_parameters( - space_station_alt_m=self.altitude, - earth_station_alt_m=self.earth_station_alt_m, - earth_station_lat_deg=self.earth_station_lat_deg, - earth_station_long_diff_deg=self.earth_station_long_diff_deg, - season=self.season - ) - if __name__ == "__main__": diff --git a/sharc/parameters/parameters_orbit.py b/sharc/parameters/parameters_orbit.py index cad76e7ee..2497c63a7 100644 --- a/sharc/parameters/parameters_orbit.py +++ b/sharc/parameters/parameters_orbit.py @@ -1,5 +1,6 @@ from dataclasses import dataclass from sharc.parameters.parameters_base import ParametersBase +from warnings import warn @dataclass @@ -17,6 +18,19 @@ class ParametersOrbit(ParametersBase): apogee_alt_km: float = 1414.0 initial_mean_anomaly: float = 0.0 + # by default, three time dependant angles are set to be + # independent random variables. + # You can model time as the only random variable instead + # by setting this parameter to True + model_time_as_random_variable: bool = False + + # You may specify a time range lower bound [s] + t_min: float = 0.0 + + # You may specify a time range upper bound [s] + # or let a default be chosen where needed + t_max: float = None + def load_parameters_from_file(self, config_file: str): """Load parameters from file and validate.""" super().load_parameters_from_file(config_file) @@ -56,3 +70,27 @@ def validate(self, ctx: str): f"Invalid {ctx}.phasing_deg = { self.phasing_deg}. \ Must be in the range [0, 360] degrees.") + + if self.enable_time_as_only_random_variable: + if self.t_max is None: + warn( + f"{ctx}.t_max was not set. Default values will be used when needed" + ) + else: + if self.t_max < self.t_min: + raise ValueError(f"{ctx}.t_max should be >= {ctx}.t_min") + if self.t_min < 0: + raise ValueError(f"{ctx}.t_min should be >= 0") + + if not self.enable_time_as_only_random_variable: + # check that defaults haven't been changed + if self.t_max != ParametersOrbit.t_max: + raise ValueError( + f"You should only set {ctx}.t_max " + f"if {ctx}.enable_time_as_only_random_variable is set to True" + ) + if self.t_min != ParametersOrbit.t_min: + raise ValueError( + f"You should only set {ctx}.t_min " + f"if {ctx}.enable_time_as_only_random_variable is set to True" + ) diff --git a/sharc/parameters/parameters_p452.py b/sharc/parameters/parameters_p452.py index 4c4179c02..c0eba01aa 100644 --- a/sharc/parameters/parameters_p452.py +++ b/sharc/parameters/parameters_p452.py @@ -37,6 +37,8 @@ class ParametersP452(ParametersBase): # determine whether clutter loss following ITU-R P.2108 is added # (TRUE/FALSE) clutter_loss: bool = True + # Determine if clutter is applied to "one_end" or "both_ends" + clutter_type: str = "one_end" def load_from_paramters(self, param: ParametersBase): """Used to load parameters of P.452 from IMT or system parameters @@ -59,3 +61,4 @@ def load_from_paramters(self, param: ParametersBase): self.rx_lat = param.rx_lat self.polarization = param.polarization self.clutter_loss = param.clutter_loss + self.clutter_type = param.clutter_type diff --git a/sharc/parameters/parameters_p619.py b/sharc/parameters/parameters_p619.py index 72d858248..81d3354c2 100644 --- a/sharc/parameters/parameters_p619.py +++ b/sharc/parameters/parameters_p619.py @@ -22,13 +22,14 @@ class ParametersP619(ParametersBase): # earth_station_long_diff_deg - difference between longitudes of IMT and satellite system # (positive if space-station is to the East of earth-station) # season - season of the year. - space_station_alt_m: float = 20000.0 earth_station_alt_m: float = 0.0 earth_station_lat_deg: float = 0.0 - earth_station_long_diff_deg: float = 0.0 season: str = "SUMMER" shadowing: bool = True noise_temperature: float = 290.0 + # Average height of clutter. According to 2018 it can be "Low", "Mid" or "High" + mean_clutter_height: str = "high" + below_rooftop: float = 100 def load_from_paramters(self, param: ParametersBase): """Used to load parameters of P.619 from IMT or system parameters @@ -38,28 +39,27 @@ def load_from_paramters(self, param: ParametersBase): param : ParametersBase IMT or system parameters """ - self.space_station_alt_m = param.space_station_alt_m self.earth_station_alt_m = param.earth_station_alt_m self.earth_station_lat_deg = param.earth_station_lat_deg - self.earth_station_long_diff_deg = param.earth_station_long_diff_deg self.season = param.season + self.mean_clutter_height = param.param_p619.mean_clutter_height if self.season.upper() not in ["SUMMER", "WINTER"]: raise ValueError(f"{self.__class__.__name__}: \ Invalid value for parameter season - {self.season}. \ Possible values are \"SUMMER\", \"WINTER\".") + allowed = {"low", "mid", "high"} + if str(self.mean_clutter_height).lower() not in allowed: + raise ValueError("Invalid type of mean_clutter_height. mean_clutter_height must be 'Low', 'Mid', or 'High'") + def set_external_parameters(self, *, - space_station_alt_m: float, earth_station_alt_m: float, earth_station_lat_deg: float, - earth_station_long_diff_deg: float, season: typing.Literal["SUMMER", "WINTER"]): """ Set external parameters for P619 propagation calculations. """ - self.space_station_alt_m = space_station_alt_m self.earth_station_alt_m = earth_station_alt_m self.earth_station_lat_deg = earth_station_lat_deg - self.earth_station_long_diff_deg = earth_station_long_diff_deg self.season = season diff --git a/sharc/parameters/parameters_single_space_station.py b/sharc/parameters/parameters_single_space_station.py index 9ddb3c749..267f87dd3 100644 --- a/sharc/parameters/parameters_single_space_station.py +++ b/sharc/parameters/parameters_single_space_station.py @@ -174,37 +174,10 @@ def propagate_parameters(self): """ Propagate relevant parameters to nested P619 and antenna objects. """ - if self.channel_model == "P619": - if self.param_p619.earth_station_alt_m != ParametersP619.earth_station_alt_m: - raise ValueError( - f"{self.section_name}.param_p619.earth_station_alt_m should not be set by hand." - "It is automatically set by other parameters in system" - ) - if self.param_p619.earth_station_lat_deg != ParametersP619.earth_station_lat_deg: - raise ValueError( - f"{self.section_name}.param_p619.earth_station_lat_deg should not be set by hand." - "It is automatically set by other parameters in system" - ) - if self.param_p619.space_station_alt_m != ParametersP619.space_station_alt_m: - raise ValueError( - f"{self.section_name}.param_p619.space_station_alt_m should not be set by hand." - "It is automatically set by other parameters in system" - ) - if self.param_p619.earth_station_lat_deg != ParametersP619.earth_station_lat_deg: - raise ValueError( - f"{self.section_name}.param_p619.earth_station_lat_deg should not be set by hand." - "It is automatically set by other parameters in system" - ) - self.param_p619.space_station_alt_m = self.geometry.altitude + self.param_p619.earth_station_alt_m = self.geometry.es_altitude self.param_p619.earth_station_lat_deg = self.geometry.es_lat_deg - if self.geometry.location.type == "FIXED": - self.param_p619.earth_station_long_diff_deg = self.geometry.location.fixed.long_deg - \ - self.geometry.es_long_deg - else: - self.param_p619.earth_station_long_diff_deg = None - # this is needed because nested parameters # don't know/cannot access parents attributes self.antenna.set_external_parameters( diff --git a/sharc/parameters/parameters_space_station.py b/sharc/parameters/parameters_space_station.py index 67df65f2d..eb7380f56 100644 --- a/sharc/parameters/parameters_space_station.py +++ b/sharc/parameters/parameters_space_station.py @@ -55,14 +55,6 @@ class ParametersSpaceStation(ParametersBase): # (positive if space-station is to the East of earth-station) # season - season of the year. param_p619 = ParametersP619() - earth_station_alt_m: float = 0.0 - earth_station_lat_deg: float = 0.0 - earth_station_long_diff_deg: float = 0.0 - season: str = "SUMMER" - - # This parameter is also used by P619, but should not be set manually. - # this should always == params.altitude - space_station_alt_m: float = None def load_parameters_from_file(self, config_file: str): """ @@ -80,11 +72,6 @@ def load_parameters_from_file(self, config_file: str): """ super().load_parameters_from_file(config_file) - if self.space_station_alt_m is not None: - raise ValueError( - "'space_station_alt_m' should not be set manually. It is always equal to 'altitude'", - ) - if self.nadir_angle != 0 and self.elevation != 0: raise ValueError( "'elevation' and 'nadir_angle' should not both be set at the same time. Choose either\ @@ -123,7 +110,6 @@ def set_derived_parameters(self): This method sets the derived parameters such as space station altitude and nadir angle based on the current configuration. It also updates the P619 parameters if applicable. """ - self.space_station_alt_m = self.altitude if self.param_p619: self.param_p619.load_from_paramters(self) diff --git a/sharc/propagation/Dataset/generate_loss_table.py b/sharc/propagation/Dataset/generate_loss_table.py index 2fac68e3c..4002987ed 100644 --- a/sharc/propagation/Dataset/generate_loss_table.py +++ b/sharc/propagation/Dataset/generate_loss_table.py @@ -11,10 +11,8 @@ # Constants frequency_MHz = 2680.0 -space_station_alt_m = 35786.0 * 1000 # Geostationary orbit altitude in meters earth_station_alt_m = 1000.0 earth_station_lat_deg = -15.7801 -earth_station_long_diff_deg = 0.0 # Not used in the loss calculation season = "SUMMER" apparent_elevation = np.arange(0, 90) # Elevation angles from 0 to 90 degrees city_name = "BRASILIA" @@ -23,10 +21,8 @@ random_number_gen = np.random.RandomState(101) propagation = PropagationP619( random_number_gen=random_number_gen, - space_station_alt_m=space_station_alt_m, earth_station_alt_m=earth_station_alt_m, earth_station_lat_deg=earth_station_lat_deg, - earth_station_long_diff_deg=earth_station_long_diff_deg, season=season, ) diff --git a/sharc/propagation/Dataset/localidades.csv b/sharc/propagation/Dataset/locations.csv similarity index 69% rename from sharc/propagation/Dataset/localidades.csv rename to sharc/propagation/Dataset/locations.csv index 151eaae6f..17be18aef 100644 --- a/sharc/propagation/Dataset/localidades.csv +++ b/sharc/propagation/Dataset/locations.csv @@ -1,2 +1,3 @@ Cidade,Latitude BRASILIA,-15.7801 +FOZ,-25.5549751 \ No newline at end of file diff --git a/sharc/propagation/propagation_clear_air_452.py b/sharc/propagation/propagation_clear_air_452.py index 187ad2180..26095ca24 100644 --- a/sharc/propagation/propagation_clear_air_452.py +++ b/sharc/propagation/propagation_clear_air_452.py @@ -1,3 +1,4 @@ + # -*- coding: utf-8 -*- """Implements ITU-R P.452 clear-air propagation model and related calculations.""" import numpy as np @@ -1824,6 +1825,7 @@ def get_loss( frequency=frequency * 1000, distance=distance * 1000, station_type=StationType.FSS_ES, + clutter_type=self.model_params.clutter_type ) else: clutter_loss = np.zeros(distance.shape) diff --git a/sharc/propagation/propagation_clutter_loss.py b/sharc/propagation/propagation_clutter_loss.py index b4d6cac6a..770bf6213 100644 --- a/sharc/propagation/propagation_clutter_loss.py +++ b/sharc/propagation/propagation_clutter_loss.py @@ -5,7 +5,6 @@ import numpy as np import scipy -import math import matplotlib.pyplot as plt @@ -75,7 +74,6 @@ def get_loss(self, *args, **kwargs) -> np.array: f = kwargs["frequency"] loc_per = kwargs.pop("loc_percentage", "RANDOM") type = kwargs["station_type"] - d = kwargs["distance"] if f.size == 1: @@ -89,55 +87,143 @@ def get_loss(self, *args, **kwargs) -> np.array: p2 = loc_per * np.ones(d.shape) if type is StationType.IMT_BS or type is StationType.IMT_UE or type is StationType.FSS_ES: - loss = self.get_terrestrial_clutter_loss( - f, d, p1, True) + self.get_terrestrial_clutter_loss(f, d, p2, False) + clutter_type = kwargs["clutter_type"] + if clutter_type == 'one_end': + loss = self.get_terrestrial_clutter_loss(f, d, p1, True) + elif clutter_type == 'both_ends': + loss = self.get_terrestrial_clutter_loss(f, d, p1, True) + self.get_terrestrial_clutter_loss(f, d, p2, False) + else: + raise ValueError("Invalid type of Clutter-type. It can be either 'one_end' or 'both-ends'") else: theta = kwargs["elevation"] - loss = self.get_spacial_clutter_loss(f, theta, p1) + earth_station_height = kwargs["earth_station_height"] + mean_clutter_height = kwargs["mean_clutter_height"] + below_rooftop = kwargs["below_rooftop"] + loss = self.get_spacial_clutter_loss(f, theta, p1, earth_station_height, mean_clutter_height) + mult_1 = np.zeros(d.shape) + num_ones = int(np.round(mult_1.size * below_rooftop / 100)) + indices = np.random.choice(mult_1.size, size=num_ones, replace=False) + mult_1.flat[indices] = 1 + loss *= mult_1 return loss def get_spacial_clutter_loss( self, frequency: float, elevation_angle: float, loc_percentage, + earth_station_height, + mean_clutter_height ): """ - This method models the calculation of the statistical distribution of - clutter loss where one end of the interference path is within man-made - clutter, and the other is a satellite, aeroplane, or other platform - above the surface of the Earth. An additional "loss" is calculated - which can be added to the basic transmission loss of a path calculated. - This model is applicable to urban and suburban environments. The method - used to develop this model is described in Report ITU-R P.2402-0. The - clutter loss not exceeded for loc_percentage% of locations "loss" for - the terrestrial to airborne or satellite path is given by this method. + Computes clutter loss according to ITU-R P.2108-2 §3.3 for Earth-space and aeronautical paths. - Parameters - ---- - frequency : center frequency [MHz] - Frequency range: 10 to 100 GHz - elevation_angle : elevation angle [degrees] - Elevation angle - range: 0 to 90 degrees - loc_percentage : percentage of locations - Percentage locations - range: 0 < p < 100 + Parameters: + - frequency: Frequency (GHz), 0.5 <= frequency <= 100 + - elevation_angle: Elevation angle (degrees), 0 <= elevation_angle <= 90 + - loc_percentage: Percentage of locations (%), 0 < loc_percentage < 100 + - earth_station_height: Ground station height (m), >= 1 + - mean_clutter_height: Median clutter height (m): Low-rise: <= 8; Mid-rise: 8 < ... <= 20; High-rise: > 20 - Returns - ------- - loss : The clutter loss not exceeded for loc_percentage% of - locations for the terrestrial to terrestrial path + Returns: + - Lces: Clutter loss according to P.2108 §3.3 (dB) """ - k1 = 93 * (frequency * 1e-3)**0.175 - A1 = 0.05 - - y = np.tan(A1 * (1 - (elevation_angle / 90)) + - math.pi * (elevation_angle / 180)) - cot = (1 / y) - - invQ = np.sqrt(2) * scipy.special.erfcinv(2 * loc_percentage) - - loss = (-k1 * (np.log(1 - loc_percentage)) * cot)**(0.5 * \ - (90 - elevation_angle) / 90) - 1 - 0.6 * invQ - - return loss + # Converting to GHz + frequency = frequency / 1000 + # Convert to percentage + loc_percentage = loc_percentage * 100 + ## Check mean_clutter_height + allowed = {"low", "mid", "high"} + if str(mean_clutter_height).lower() not in allowed: + raise ValueError("Invalid type of mean_clutter_height. mean_clutter_height must be 'Low', 'Mid', or 'High'") + # --- Table 7: plos parameters --- + if str(mean_clutter_height).lower() == "low": + ak, bk, ck = 4.9, 6.7, 2.6 + aC, bC = 0.19, 0 + aV, bV = 1.4, 74 + elif str(mean_clutter_height).lower() == "mid": + ak, bk, ck = -2.6, 6.6, 2.0 + aC, bC = 0.42, -6.7 + aV, bV = 0.15, 97 + else: + ak, bk, ck = 2.4, 7.0, 1.0 + aC, bC = 0.19, -2.7 + aV, bV = 0.15, 98 + + # --- Table 8: p_fclo_los parameters --- + if str(mean_clutter_height).lower() == "low": + akp, bkp = 6.0, 0.07 + aCp, bC1p, bC2p = 0.15, 5.4, -0.3 + bC3p, bC4p = 3.2, 0.07 + cCp, aVp, bVp = -27, 1.6, -17 + elif str(mean_clutter_height).lower() == "mid": + akp, bkp = 3.6, 0.05 + aCp, bC1p, bC2p = 0.17, 13, -0.2 + bC3p, bC4p = 3.7, 0.05 + cCp, aVp, bVp = -41, 1, -21 + else: + akp, bkp = 5.0, 0.003 + aCp, bC1p, bC2p = 0.17, 32.6, 0.012 + bC3p, bC4p = -23.9, -0.07 + cCp, aVp, bVp = -41, 1, -18 + + # --- LoS probability (eqs. 7-8) --- + Vmax = np.minimum(aV * earth_station_height + bV, 100) + Ce = aC * earth_station_height + bC + k = ((earth_station_height + ak) / bk) ** ck + num = 1 - np.exp(-k * (elevation_angle + Ce) / 90) + den = 1 - np.exp(-k * (90 + Ce) / 90) + pLoS = np.maximum(0, Vmax * (num / den)) + + # --- Conditional probability of Fresnel clearance (eqs. 9-10) --- + Vmaxp = np.minimum(aVp * earth_station_height + bVp, 0) * frequency ** (-0.55) + 100 + Cep = (frequency * 1e9) ** aCp + bC1p * np.exp(bC2p * earth_station_height) + bC3p * np.exp(bC4p * + earth_station_height) + cCp + kp = akp * np.exp(bkp * earth_station_height) + nump = 1 - np.exp(-kp * (elevation_angle + Cep) / 90) + denp = 1 - np.exp(-kp * (90 + Cep) / 90) + pFcLoS_LoS = np.maximum(0, Vmaxp * (nump / denp)) + pFcLoS = pLoS * pFcLoS_LoS / 100 + + # Table 9 constants + alpha1, alpha2 = 8.54, 0.056 + beta1, beta2 = 17.57, 6.32 + gamma1, gamma2 = 0.63, 0.19 + + mu = alpha1 + beta1 * np.log(1 + (90 - elevation_angle) / 90) + frequency ** gamma1 + sigma = alpha2 + beta2 * np.log(1 + (90 - elevation_angle) / 90) + frequency ** gamma2 + + # Table 10 constants + theta1, theta2 = -3.542, -155.1 + sigma1, sigma2 = -1.06, -6.342 + + # Prepare output array + Lces = np.zeros_like(frequency, dtype=float) + + # Branch logic: mask arrays for each regime + mask_nlos = loc_percentage > pLoS + mask_fclo = loc_percentage <= pFcLoS + mask_between = (~mask_nlos) & (~mask_fclo) + + # NLoS branch + if np.any(mask_nlos): + pp = (loc_percentage[mask_nlos] - pLoS[mask_nlos]) / (100 - pLoS[mask_nlos]) + Finv = -np.sqrt(2) * scipy.special.erfcinv(2 * pp) + Lces_nlos = mu[mask_nlos] + sigma[mask_nlos] * Finv + Lces[mask_nlos] = np.maximum(Lces_nlos, 6) + + # Fresnel-clear branch + if np.any(mask_fclo): + ratio = loc_percentage[mask_fclo] / pFcLoS[mask_fclo] + Lces[mask_fclo] = (theta1 * np.exp(theta2 * ratio) + sigma1 * np.exp(sigma2 * ratio)) + + # Between + if np.any(mask_between): + Lces[mask_between] = ( + 6 * (loc_percentage[mask_between] - pFcLoS[mask_between]) / + (pLoS[mask_between] - pFcLoS[mask_between]) + ) + + return Lces def get_terrestrial_clutter_loss( self, @@ -192,12 +278,10 @@ def get_terrestrial_clutter_loss( loss_2km = np.zeros(d.shape) if apply_both_ends: - # minimum path length for the correction to be applied at only one - # end of the path + # minimum path length for the correction to be applied at only one end of the path id_d = np.where(d >= 1000)[0] else: - # minimum path length for the correction to be applied at both ends - # of the path + # minimum path length for the correction to be applied at both ends of the path id_d = np.where(d >= 250)[0] if len(id_d): @@ -207,17 +291,16 @@ def get_terrestrial_clutter_loss( Ls_temp = 32.98 + 3.0 * np.log10(f[id_d] * 1e-3) Ls = 23.9 * np.log10(d[id_d] * 1e-3) + Ls_temp invQ = np.sqrt(2) * scipy.special.erfcinv(2 * (p[id_d])) - sigma_cb = np.sqrt(((sigma_l**(2.0)) * (10.0**(-0.2 * Ll)) + (sigma_s**( - 2.0)) * (10.0**(-0.2 * Ls))) / (10.0**(-0.2 * Ll) + 10.0**(-0.2 * Ls))) + sigma_cb = np.sqrt(((sigma_l**(2.0)) * (10.0**(-0.2 * Ll)) + (sigma_s**(2.0)) * + (10.0**(-0.2 * Ls))) / (10.0**(-0.2 * Ll) + 10.0**(-0.2 * Ls))) loss[id_d] = -5.0 * \ np.log10(10 ** (-0.2 * Ll) + 10 ** (-0.2 * Ls)) - sigma_cb * invQ - # The clutter loss must not exceed a maximum value calculated for 𝑑 - # = 2 𝑘m (loss_2km) + # The clutter loss must not exceed a maximum value calculated for 𝑑 = 2 𝑘m (loss_2km) Ls_2km = 23.9 * np.log10(2) + Ls_temp - sigma_cb_2km = np.sqrt(((sigma_l**(2.0)) * (10.0**(-0.2 * Ll)) + (sigma_s**( - 2.0)) * (10.0**(-0.2 * Ls_2km))) / (10.0**(-0.2 * Ll) + 10.0**(-0.2 * Ls_2km))) + sigma_cb_2km = np.sqrt(((sigma_l**(2.0)) * (10.0**(-0.2 * Ll)) + (sigma_s**(2.0)) * + (10.0**(-0.2 * Ls_2km))) / (10.0**(-0.2 * Ll) + 10.0**(-0.2 * Ls_2km))) loss_2km[id_d] = -5.0 * \ np.log10(10 ** (-0.2 * Ll) + 10 ** (-0.2 * Ls_2km)) - \ sigma_cb_2km * invQ @@ -232,32 +315,34 @@ def get_terrestrial_clutter_loss( if __name__ == '__main__': elevation_angle = np.array([90, 80, 70, 60, 50, 40, 30, 20, 15, 10, 5, 0]) - loc_percentage = np.linspace(0.01, 0.995, 1000) - frequency = 30000 * np.ones(elevation_angle.shape) - + loc_percentage = np.linspace(0.1 / 100, 99.9 / 100, 1001) + frequency = 3000 * np.ones(loc_percentage.shape) + earth_station_height = 5 * np.ones(loc_percentage.shape) random_number_gen = np.random.RandomState(101) cl = PropagationClutterLoss(random_number_gen) clutter_loss = np.empty([len(elevation_angle), len(loc_percentage)]) - for i in range(len(loc_percentage)): - clutter_loss[:, i] = cl.get_spacial_clutter_loss( + for i in range(len(elevation_angle)): + clutter_loss[i, :] = cl.get_spacial_clutter_loss( frequency, - elevation_angle, - loc_percentage[i], + elevation_angle[i], + loc_percentage, + earth_station_height, + 'Low', ) fig = plt.figure(figsize=(8, 6), facecolor='w', edgecolor='k') ax = fig.gca() for j in range(len(elevation_angle)): - ax.plot(clutter_loss[j, :], 100 * loc_percentage, + ax.plot(clutter_loss[j, :], loc_percentage, label="%i deg" % elevation_angle[j], linewidth=1) plt.title("Cumulative distribution of clutter loss not exceeded for 30 GHz") plt.xlabel("clutter loss [dB]") plt.ylabel("percent of locations [%]") - plt.xlim((-10, 70)) + plt.xlim((-5, 30)) plt.ylim((0, 100)) plt.legend(loc="lower right") plt.tight_layout() diff --git a/sharc/propagation/propagation_factory.py b/sharc/propagation/propagation_factory.py index a40679d30..4f99144ea 100644 --- a/sharc/propagation/propagation_factory.py +++ b/sharc/propagation/propagation_factory.py @@ -79,7 +79,7 @@ def create_propagation( param.imt.topology.type}", ) else: # P.619 model is used only for space-to-earth links - if param.imt.topology.type != "NTN" and not param_system.is_space_to_earth: + if param.imt.topology.type not in ["NTN", "MSS_DC"] and not param_system.is_space_to_earth: raise ValueError(( "PropagationFactory: Channel model P.619 is invalid" f"for system {param.general.system} and IMT " @@ -87,12 +87,13 @@ def create_propagation( )) return PropagationP619( random_number_gen=random_number_gen, - space_station_alt_m=param_system.param_p619.space_station_alt_m, earth_station_alt_m=param_system.param_p619.earth_station_alt_m, earth_station_lat_deg=param_system.param_p619.earth_station_lat_deg, - earth_station_long_diff_deg=param_system.param_p619.earth_station_lat_deg, - season=param_system.season, + season=param_system.param_p619.season, + mean_clutter_height=param_system.param_p619.mean_clutter_height, + below_rooftop=param_system.param_p619.below_rooftop ) + elif channel_model == "P452": return PropagationClearAir( random_number_gen, param_system.param_p452) diff --git a/sharc/propagation/propagation_p619.py b/sharc/propagation/propagation_p619.py index 2ef7f81c7..a5c4907da 100644 --- a/sharc/propagation/propagation_p619.py +++ b/sharc/propagation/propagation_p619.py @@ -11,6 +11,7 @@ from scipy.interpolate import interp1d import numpy as np from multipledispatch import dispatch +from warnings import warn from sharc.station_manager import StationManager from sharc.parameters.parameters import Parameters from sharc.propagation.propagation import Propagation @@ -34,11 +35,11 @@ class PropagationP619(Propagation): def __init__( self, random_number_gen: np.random.RandomState, - space_station_alt_m: float, earth_station_alt_m: float, earth_station_lat_deg: float, - earth_station_long_diff_deg: float, season: str, + mean_clutter_height: str, + below_rooftop: float ): """Implements the earth-to-space channel model from ITU-R P.619 @@ -46,8 +47,6 @@ def __init__( ---------- random_number_gen : np.random.RandomState randon number generator - space_station_alt_m : float - The space station altitude in meters earth_station_alt_m : float The Earth station altitude in meters earth_station_lat_deg : float @@ -77,11 +76,11 @@ def __init__( self.surf_water_dens_has_atmospheric_loss = [] self.atmospheric_loss = [] self.elevation_delta = .01 + self.mean_clutter_height = mean_clutter_height + self.below_rooftop = below_rooftop - self.space_station_alt_m = space_station_alt_m self.earth_station_alt_m = earth_station_alt_m self.earth_station_lat_deg = earth_station_lat_deg - self.earth_station_long_diff_deg = earth_station_long_diff_deg if season.upper() not in ["SUMMER", "WINTER"]: raise ValueError( @@ -93,10 +92,12 @@ def __init__( self.city_name = self._get_city_name_by_latitude() if self.city_name != "Unknown": self.lookup_table = True + else: + warn('Using analytical model for atmospheric attenuation. No lookup table available for this latitude.') def _get_city_name_by_latitude(self): localidades_file = os.path.join( - os.path.dirname(__file__), 'Dataset/localidades.csv', + os.path.dirname(__file__), 'Dataset/locations.csv', ) with open(localidades_file, mode='r') as file: reader = csv.DictReader(file) @@ -130,12 +131,9 @@ def _get_atmospheric_gasses_loss(self, *args, **kwargs) -> float: if lookupTable and self.city_name != 'Unknown': # Define the path to the CSV file output_dir = os.path.join(os.path.dirname(__file__), 'Dataset') + lookup_table_name = f'{self.city_name}_{int(frequency_MHz)}_{int(self.earth_station_alt_m)}m.csv' csv_file = os.path.join( - output_dir, f'{ - self.city_name}_{ - int(frequency_MHz)}_{ - int( - self.earth_station_alt_m)}m.csv', ) + output_dir, lookup_table_name) if os.path.exists(csv_file): elevations = [] losses = [] @@ -148,6 +146,10 @@ def _get_atmospheric_gasses_loss(self, *args, **kwargs) -> float: interpolation_function = interp1d( elevations, losses, kind='linear', fill_value='extrapolate', ) return interpolation_function(apparent_elevation) + else: + raise FileNotFoundError( + f"CSV file {lookup_table_name} not found but lookupTable is set to True. Did you configured the 'Dataset/locations.csv' lookup table correctly? ", + ) earth_radius_km = EARTH_RADIUS / 1000 a_acc = 0. # accumulated attenuation (in dB) @@ -348,6 +350,7 @@ def get_loss( # Elevation angles seen from the station on Earth. elevation_angles = {} if station_a.is_space_station: + earth_station_height = station_b.height elevation_angles["free_space"] = station_b.get_elevation(station_a) earth_station_antenna_gain = station_b_gains # if (station_b_gains.shape != distance.shape): @@ -362,6 +365,7 @@ def get_loss( elevation_angles["apparent"] = np.transpose( elevation_angles["apparent"]) elif station_b.is_space_station: + earth_station_height = station_a.height elevation_angles["free_space"] = station_a.get_elevation(station_b) earth_station_antenna_gain = station_a_gains elevation_angles["apparent"] = self.apparent_elevation_angle( @@ -408,11 +412,12 @@ def get_loss( is_earth_to_space_link, earth_station_antenna_gain, is_single_entry_interf, + earth_station_height, ) return loss - @dispatch(np.ndarray, np.ndarray, np.ndarray, dict, bool, np.ndarray, bool) + @dispatch(np.ndarray, np.ndarray, np.ndarray, dict, bool, np.ndarray, bool, np.ndarray) def get_loss( self, distance: np.array, @@ -422,6 +427,8 @@ def get_loss( earth_to_space: bool, earth_station_antenna_gain: np.array, single_entry: bool, + earth_station_height: np.array, + ) -> np.array: """ Calculates path loss for earth-space link @@ -481,6 +488,9 @@ def get_loss( distance=distance, elevation=elevation["free_space"], station_type=StationType.FSS_SS, + earth_station_height=earth_station_height, + mean_clutter_height=self.mean_clutter_height, + below_rooftop=self.below_rooftop, ) building_loss = self.building_entry.get_loss( frequency, elevation["apparent"], @@ -505,31 +515,17 @@ def get_loss( import matplotlib.pyplot as plt - # params = Parameters() - - # propagation_path = os.getcwd() - # sharc_path = os.path.dirname(propagation_path) - # param_file = os.path.join(sharc_path, "parameters", "parameters.ini") - # deprecated - - # params.set_file_name(param_file) - # params.read_params() - - # sat_params = params.fss_ss - - space_station_alt_m = 20000.0 earth_station_alt_m = 1000.0 earth_station_lat_deg = -15.7801 - earth_station_long_diff_deg = 0.0 season = "SUMMER" random_number_gen = np.random.RandomState(101) propagation = PropagationP619( random_number_gen=random_number_gen, - space_station_alt_m=space_station_alt_m, earth_station_alt_m=earth_station_alt_m, earth_station_lat_deg=earth_station_lat_deg, - earth_station_long_diff_deg=earth_station_long_diff_deg, + mean_clutter_height='low', + below_rooftop=0.0, season=season, ) diff --git a/sharc/satellite/ngso/custom_functions.py b/sharc/satellite/ngso/custom_functions.py index cd20e65f4..6442a35c1 100644 --- a/sharc/satellite/ngso/custom_functions.py +++ b/sharc/satellite/ngso/custom_functions.py @@ -68,7 +68,7 @@ def keplerian2eci(a, e, delta, Omega, omega, nu): delta : float Orbital inclination in degrees. Omega : np.ndarray - Right ascension of the ascending node (RAAN) in degrees. Shape = (N,) + Right ascension of the ascending node (RAAN) in degrees. Shape = (N, length(t)) omega : float Argument of perigee in degrees. nu : np.ndarray @@ -92,10 +92,8 @@ def keplerian2eci(a, e, delta, Omega, omega, nu): # Trigonometric calculations cos_gamma = np.cos(gamma) sin_gamma = np.sin(gamma) - # Shape (N, 1) for broadcasting - cos_raan = np.cos(-Omega_rad)[:, np.newaxis] - # Shape (N, 1) for broadcasting - sin_raan = np.sin(-Omega_rad)[:, np.newaxis] + cos_raan = np.cos(-Omega_rad) + sin_raan = np.sin(-Omega_rad) cos_incl = np.cos(delta_rad) sin_incl = np.sin(delta_rad) diff --git a/sharc/satellite/ngso/orbit_model.py b/sharc/satellite/ngso/orbit_model.py index d97f13660..92f598394 100644 --- a/sharc/satellite/ngso/orbit_model.py +++ b/sharc/satellite/ngso/orbit_model.py @@ -19,7 +19,12 @@ def __init__(self, delta: float, hp: float, ha: float, - Mo: float): + Mo: float, + *, + model_time_as_random_variable: bool, + t_min: float, + t_max: float | None + ): """Instantiates and OrbitModel object from the Orbit parameters as specified in S.1529. Parameters @@ -42,11 +47,20 @@ def __init__(self, altitude of apogee in km Mo : float initial mean anomaly for first satellite of first plane, in degrees + model_time_as_random_variable: bool + whether get_orbit_positions_random will use only time as random variable + t_min: float + if model_time_as_random_variable == True, + defines the lower bound of the time distribution + t_max: float + if model_time_as_random_variable == True, + defines the upper bound of the time distribution """ self.Nsp = Nsp self.Np = Np self.phasing = phasing self.long_asc = long_asc + self.omega_0 = omega self.omega = omega self.delta = delta self.perigee_alt_km = hp @@ -58,12 +72,14 @@ def __init__(self, hp + ha + 2 * EARTH_RADIUS_KM) / 2 # semi-major axis, in km self.eccentricity = (ha + EARTH_RADIUS_KM - self.semi_major_axis) / \ self.semi_major_axis # orbital eccentricity (e) - self.orbital_period_sec = 2 * np.pi * \ - np.sqrt(self.semi_major_axis ** 3 / KEPLER_CONST) # orbital period, in seconds + # TODO: consider J2 perturbations for mean motion and orbital period + self.mean_motion = np.sqrt(KEPLER_CONST / self.semi_major_axis ** 3) + # orbital period, in seconds + self.orbital_period_sec = 2 * np.pi / self.mean_motion # satellite separation angle in the plane (degrees) self.sat_sep_angle_deg = 360 / self.Nsp # angle between plane intersections with the equatorial plane (degrees) - self.orbital_plane_inclination = 360 / self.Np + self.orbital_plane_spacing = 360 / self.Np # Initial mean anomalies for all the satellites initial_mean_anomalies_deg = ( @@ -76,9 +92,21 @@ def __init__(self, # Initial longitudes of ascending node for all the planes inital_raan = (self.long_asc * np.ones(self.Nsp) + np.arange(self.Np) - [:, np.newaxis] * self.orbital_plane_inclination) % 360 + [:, np.newaxis] * self.orbital_plane_spacing) % 360 self.inital_raan_rad = np.radians(inital_raan.flatten()) + self.model_time_as_random_variable = model_time_as_random_variable + self.t_min = t_min + + if t_max is not None: + self.t_max = t_max + else: + # sane default + self.t_max = t_min + self.orbital_period_sec * 1e3 + # TODO: implement 1503 secular drift to enable correct higher time ceiling + # self.t_max = t_min + 365 * 24 * 60 * 60 + # self.t_max = sys.float_info.max + # Calculated variables self.mean_anomaly = None # computed mean anomalies self.raan_rad = None # computed longitudes of the ascending node @@ -117,12 +145,52 @@ def get_satellite_positions_time_interval( self.orbital_period_sec + interval_secs, interval_secs) - return self.__get_satellite_positions(t) + return self.get_orbit_positions_time_instant(t) def get_orbit_positions_time_instant(self, time_instant_secs=0) -> dict: - """Returns satellite positions in a determined time instant in seconds""" - t = np.array([time_instant_secs]) - return self.__get_satellite_positions(t) + """Returns the Satellite positins for a given time vector within the orbit period. + + Parameters + ---------- + time_instant_secs: np.array + time instants inside the orbit period in seconds + + Returns + ------- + dict + A dictionary with satellite positions in spherical and ecef coordinates. + lat, lon, alt, sx, sy, sz + """ + t = np.atleast_1d(time_instant_secs) + + assert t.ndim == 1, "Input must be scalar or 1D array" + + # TODO: add J2 perturbations based on 1503 + # Mean anomaly (M) + self.mean_anomaly = ( + self.initial_mean_anomalies_rad[:, None] + self.mean_motion * t + ) % (2 * np.pi) + + # TODO: add J2 perturbations based on 1503 + # Longitudes of the ascending node (OmegaG) + # shape (Np*Nsp, len(t)) + self.raan_rad = self.inital_raan_rad[:, None] + + # TODO: add J2 perturbations based on 1503 + # perigee argument + self.omega = self.omega_0 + + self.omega_rad = np.deg2rad(self.omega) + + # The time to be used for calculation of Earth's rotation angle + earth_rotated_t = t + + return self.__get_satellite_positions_from_angles( + self.mean_anomaly, + self.raan_rad, + self.omega_rad, + earth_rotated_t, + ) def get_orbit_positions_random( self, @@ -141,98 +209,68 @@ def get_orbit_positions_random( A dictionary with satellite positions in spherical and ecef coordinates. lat, lon, sx, sy, sz """ - - # return self.__get_satellite_positions(rng.random_sample(1) * 1000 * self.orbital_period_sec) + if self.model_time_as_random_variable: + return self.get_orbit_positions_time_instant( + self.t_min + (self.t_max - self.t_min) * rng.random_sample(n_samples) + ) # Mean anomaly (M) self.mean_anomaly = (self.initial_mean_anomalies_rad[:, None] + 2 * np.pi * rng.random_sample(n_samples)) % (2 * np.pi) - # Eccentric anomaly (E) - self.eccentric_anom = eccentric_anomaly( - self.eccentricity, self.mean_anomaly) - - # True anomaly (v) - self.true_anomaly = 2 * np.arctan(np.sqrt((1 + self.eccentricity) / ( - 1 - self.eccentricity)) * np.tan(self.eccentric_anom / 2)) - - self.true_anomaly = np.mod(self.true_anomaly, 2 * np.pi) - - # Distance of the satellite to Earth's center (r) - r = self.semi_major_axis * \ - (1 - self.eccentricity ** 2) / (1 + self.eccentricity * np.cos(self.true_anomaly)) - - # True anomaly relative to the line of nodes (gamma) - self.true_anomaly_rel_line_nodes = wrap2pi( - self.true_anomaly + - np.radians( - self.omega)) # gamma in the interval [-pi, pi] - - # Latitudes of the satellites, in radians (theta) - # theta = np.arcsin(np.sin(gamma) * np.sin(np.radians(self.delta))) - - # Longitude variation due to angular displacement, in radians (phiS) - # phiS = np.arccos(np.cos(gamma) / np.cos(theta)) * np.sign(gamma) + # just selecting a random earth rotation for later coordinate transformation + earth_rotated_t = 2 * np.pi * rng.random_sample(n_samples) / EARTH_ROTATION_RATE # Longitudes of the ascending node (OmegaG) # shape (Np*Nsp, len(t)) - raan_rad = (self.inital_raan_rad[:, None] + - 2 * np.pi * rng.random_sample(n_samples)) - raan_rad = wrap2pi(raan_rad) - - # POSITION CALCULATION IN ECEF COORDINATES - ITU-R S.1503 - r_eci = keplerian2eci(self.semi_major_axis, - self.eccentricity, - self.delta, - np.degrees(self.inital_raan_rad), - self.omega, - np.degrees(self.true_anomaly)) - - r_ecef = eci2ecef( - 2 * - np.pi * - rng.random_sample(n_samples) / - EARTH_ROTATION_RATE, - r_eci) - sx, sy, sz = r_ecef[0], r_ecef[1], r_ecef[2] - lat = np.degrees(np.arcsin(sz / r)) - lon = np.degrees(np.arctan2(sy, sx)) - # (lat, lon, _) = ecef2lla(sx, sy, sz) - - pos_vector = { - 'lat': lat, - 'lon': lon, - 'alt': r - EARTH_RADIUS_KM, - 'sx': sx, - 'sy': sy, - 'sz': sz - } - return pos_vector - - def __get_satellite_positions(self, t: np.array) -> dict: - """Returns the Satellite positins (both lat long and ecef) for a given time vector within the orbit period. - - Parameters - ---------- - t : np.array - time instants inside the orbit period in seconds + # FIXME: previous implementation, due to unexpected behavior, did not + # use the drawn random samples. Choose whether to maintain behavior + # self.raan_rad = (self.inital_raan_rad[:, None] + + # 2 * np.pi * rng.random_sample(n_samples)) + self.raan_rad = self.inital_raan_rad[:, None] + + # perigee argument + # NOTE: using fixed perigee argument for circular orbits always make sense + self.omega = self.omega_0 + + self.omega_rad = np.deg2rad(self.omega) + + return self.__get_satellite_positions_from_angles( + self.mean_anomaly, + self.raan_rad, + self.omega_rad, + earth_rotated_t, + ) - Returns - ------- - dict - A dictionary with satellite positions in spherical and ecef coordinates. - lat, lon, sx, sy, sz + def __get_satellite_positions_from_angles( + self, + mean_anomaly: np.ndarray, + raan_rad: np.ndarray, + omega_rad: np.array, + earth_rotated_t: np.ndarray, + ): """ - # Mean anomaly (M) - self.mean_anomaly = (self.initial_mean_anomalies_rad[:, None] + ( - 2 * np.pi / self.orbital_period_sec) * t) % (2 * np.pi) + mean_anomaly: + Mean anomaly (M) + raan_rad: np.ndarray + Longitudes of the ascending node (OmegaG) + omega_rad: np.ndarray + Perigee argument + NOTE: doesn't matter for circular orbits + earth_rotated_t: + The time to be used for calculation of Earth's rotation angle + """ + assert (raan_rad.shape == (self.Np * self.Nsp, len(earth_rotated_t))) or ( + raan_rad.shape == (self.Np * self.Nsp, 1) + ) # Eccentric anomaly (E) self.eccentric_anom = eccentric_anomaly( - self.eccentricity, self.mean_anomaly) + self.eccentricity, mean_anomaly) # True anomaly (v) self.true_anomaly = 2 * np.arctan(np.sqrt((1 + self.eccentricity) / ( 1 - self.eccentricity)) * np.tan(self.eccentric_anom / 2)) + self.true_anomaly = np.mod(self.true_anomaly, 2 * np.pi) # Distance of the satellite to Earth's center (r) @@ -242,8 +280,7 @@ def __get_satellite_positions(self, t: np.array) -> dict: # True anomaly relative to the line of nodes (gamma) self.true_anomaly_rel_line_nodes = wrap2pi( self.true_anomaly + - np.radians( - self.omega)) # gamma in the interval [-pi, pi] + omega_rad) # gamma in the interval [-pi, pi] # Latitudes of the satellites, in radians (theta) # theta = np.arcsin(np.sin(gamma) * np.sin(np.radians(self.delta))) @@ -251,20 +288,17 @@ def __get_satellite_positions(self, t: np.array) -> dict: # Longitude variation due to angular displacement, in radians (phiS) # phiS = np.arccos(np.cos(gamma) / np.cos(theta)) * np.sign(gamma) - # Longitudes of the ascending node (OmegaG) - # shape (Np*Nsp, len(t)) - raan_rad = (self.inital_raan_rad[:, None] + EARTH_ROTATION_RATE * t) raan_rad = wrap2pi(raan_rad) # POSITION CALCULATION IN ECEF COORDINATES - ITU-R S.1503 r_eci = keplerian2eci(self.semi_major_axis, self.eccentricity, self.delta, - np.degrees(self.inital_raan_rad), - self.omega, + np.degrees(raan_rad), + np.degrees(omega_rad), np.degrees(self.true_anomaly)) - r_ecef = eci2ecef(t, r_eci) + r_ecef = eci2ecef(earth_rotated_t, r_eci) sx, sy, sz = r_ecef[0], r_ecef[1], r_ecef[2] lat = np.degrees(np.arcsin(sz / r)) lon = np.degrees(np.arctan2(sy, sx)) @@ -273,6 +307,7 @@ def __get_satellite_positions(self, t: np.array) -> dict: pos_vector = { 'lat': lat, 'lon': lon, + 'alt': r - EARTH_RADIUS_KM, 'sx': sx, 'sy': sy, 'sz': sz @@ -298,7 +333,10 @@ def main(): "delta": 53, "hp": 525, "ha": 525, - "Mo": 0 + "Mo": 0, + "model_time_as_random_variable": False, + "t_min": 0.0, + "t_max": None } print("Orbit parameters:") diff --git a/sharc/satellite/scripts/plot_satellites_3d.py b/sharc/satellite/scripts/plot_satellites_3d.py index a0f82748c..04f943c8c 100644 --- a/sharc/satellite/scripts/plot_satellites_3d.py +++ b/sharc/satellite/scripts/plot_satellites_3d.py @@ -20,7 +20,10 @@ delta=52, # inclination in degrees hp=1414, # perigee altitude in km ha=1414, # apogee altitude in km - Mo=0 # mean anomaly in degrees + Mo=0, # mean anomaly in degrees + model_time_as_random_variable=False, + t_min=0.0, + t_max=None ) # Set parameters for ground station diff --git a/sharc/satellite/scripts/plot_time_domain_sampling.py b/sharc/satellite/scripts/plot_time_domain_sampling.py new file mode 100644 index 000000000..2e4899b2b --- /dev/null +++ b/sharc/satellite/scripts/plot_time_domain_sampling.py @@ -0,0 +1,286 @@ +import numpy as np +import plotly.graph_objects as go + +from sharc.satellite.ngso.orbit_model import OrbitModel +from sharc.satellite.utils.sat_utils import calc_elevation, lla2ecef +from sharc.satellite.scripts.plot_globe import plot_globe_with_borders +from sharc.satellite.ngso.constants import EARTH_RADIUS_KM + + +def plot_orbit_trace(orbit): + """ + Traces the orbit satellite paths from t_min to t_max + """ + # Plot satellite traces from a time interval + fig = plot_globe_with_borders(True, None, True) + step = 5 + pos_vec = orbit.get_orbit_positions_time_instant( + np.linspace(orbit.t_min, orbit.t_max, int((orbit.t_max - orbit.t_min) / step)) + ) + fig.add_trace(go.Scatter3d(x=pos_vec['sx'].flatten(), + y=pos_vec['sy'].flatten(), + z=pos_vec['sz'].flatten(), + mode='lines', + showlegend=False)) + + fig.update_layout( + title_text='Satellite Traces', + scene=dict( + xaxis_title='X [km]', + yaxis_title='Y [km]', + zaxis_title='Z [km]', + aspectmode='data' + ), + margin=dict(l=0, r=0, b=0, t=0) + ) + + return fig + + +def plot_sampling( + orbit, + ground_sta_lat, + ground_sta_lon, + ground_sta_alt, + min_elev_angle_deg, + num_drops=1000 +): + """ + Gets num_drops random snapshots from orbit model and plots the result + """ + # Show visible satellites from ground-station + fig = plot_globe_with_borders(True, None, True) + rng = np.random.RandomState(seed=6) + pos_vec = orbit.get_orbit_positions_random(rng=rng, n_samples=num_drops) + look_angles = calc_elevation( + ground_sta_lat, + pos_vec['lat'], + ground_sta_lon, + pos_vec['lon'], + sat_height=orbit.apogee_alt_km * 1e3, + es_height=ground_sta_alt) + + # plot all satellites in drops + fig.add_trace(go.Scatter3d(x=pos_vec['sx'].flatten(), + y=pos_vec['sy'].flatten(), + z=pos_vec['sz'].flatten(), + mode='markers', + marker=dict(size=2, + color='red', + opacity=0.8), + showlegend=False)) + + # plot visible satellites + fig.add_trace(go.Scatter3d( + x=pos_vec['sx'][np.where(look_angles > min_elev_angle_deg)].flatten(), + y=pos_vec['sy'][np.where(look_angles > min_elev_angle_deg)].flatten(), + z=pos_vec['sz'][np.where(look_angles > min_elev_angle_deg)].flatten(), + mode='markers', + marker=dict(size=2, + color='green', + opacity=0.8), + showlegend=False)) + + # plot ground station + groud_sta_pos = lla2ecef(ground_sta_lat, ground_sta_lon, ground_sta_alt) + fig.add_trace(go.Scatter3d( + x=np.array(groud_sta_pos[0] / 1e3), + y=np.array(groud_sta_pos[1] / 1e3), + z=np.array(groud_sta_pos[2] / 1e3), + mode='markers', + marker=dict(size=4, + color='blue', + opacity=1.0), + showlegend=False)) + + return fig + + +def get_times_for_sat_in_correct_latitude( + earth_station_lat: np.ndarray, + orbit: OrbitModel, +): + """ + Returns time: nd.array + A time array for the times where a satellite is at the given latitude + in the [0, orbit_period_sec) range. + Since every satellite passes through that lat twice, + ascending and descending, time.shape == (Np * Nsp * 2,) + """ + inclin = np.abs(orbit.delta) + + if np.abs(earth_station_lat) > inclin: + # return None + raise ValueError( + "It is impossible to find a time where the satellite " + f"passes ahead of a point lat={earth_station_lat} if " + f"its inclination is {inclin}" + ) + # considering plane geometry, latitude of the satellite follows + # is the argument of latitude u + # u(t) = v(t) + omega(t) + # simplifying argument of perigee omega(t) = omega can be done for circular orbits + # or when ideal keplerian model is considered (since there is no precession) + # so u(t) = v(t) + omega + + # also, sin(inclin) * sin(u) = sin(lat) + # so, to get wanted time t_w, + # u = arcsin(sin(lat) / sin(inclin)) + u1 = np.arcsin( + np.sin(np.deg2rad(earth_station_lat)) / np.sin(np.deg2rad(orbit.delta)) + ) + u2 = np.pi - u1 + u = np.array([u1, u2]) + # true anomaly: + # NOTE: this would not be correct when considering J2 perigee argument precession + # on non-circular orbits. In that case, numerical methods are needed + v = u - np.deg2rad(orbit.omega_0) + + # eccentric anomaly + E = 2 * np.arctan2( + np.sqrt(1 - orbit.eccentricity) * np.sin(v / 2), + np.sqrt(1 + orbit.eccentricity) * np.cos(v / 2) + ) + + # mean anomaly wanted, to get desired latitude + M_w = E - orbit.eccentricity * np.sin(E) + + # time wanted for desired latitude is t_w (mod orbit_period) + t_w = ( + (M_w[None, :] - orbit.initial_mean_anomalies_rad[:, None]) / + orbit.mean_motion + ) % orbit.orbital_period_sec + + return t_w.flatten() + + +def min_and_max_t_from( + selected_t: float, + selected_sat_i: float, + orbit: OrbitModel, + ground_lat: float, + ground_lon: float, + ground_alt: float, + min_elevation, + step: float = 5 +): + """ + Searches and returns min and max time for which satellite is still visible. + selected_t is the midpoint of the search domain + """ + # check if it is in bounding box + low = selected_t - orbit.orbital_period_sec / 4 + high = selected_t + orbit.orbital_period_sec / 4 + ts = np.arange(low, high, step) + all_pos = orbit.get_orbit_positions_time_instant(ts) + sat_lats = all_pos["lat"][selected_sat_i] + sat_lons = all_pos["lon"][selected_sat_i] + sat_alts = all_pos["alt"][selected_sat_i] + + elev = calc_elevation( + ground_lat, sat_lats, + ground_lon, sat_lons, + sat_height=sat_alts * 1e3, + es_height=ground_alt + ) + mask = elev >= min_elevation + ts_within = np.where(mask)[0] + + t_min_i = ts_within[0] - 1 + t_max_i = ts_within[-1] + 1 + + return ts[t_min_i], ts[t_max_i] + + +if __name__ == "__main__": + orbit = OrbitModel( + Nsp=120, # number of sats per plane + Np=28, # number of planes + phasing=1.5, # phasing in degrees + long_asc=0, # longitude of the ascending node in degrees + omega=0, # argument of perigee in degrees + delta=52, # inclination in degrees + hp=525, # perigee altitude in km + ha=525, # apogee altitude in km + Mo=0, # mean anomaly in degrees + model_time_as_random_variable=True, + t_min=0.0, + t_max=None + ) + + # Set parameters for ground station + GROUND_STA_LAT = -15.7801 + GROUND_STA_LON = -42.9292 + GROUND_STA_ALT = 1200 + + MIN_ELEV_ANGLE_DEG = 5.0 + + possible_ts = get_times_for_sat_in_correct_latitude( + GROUND_STA_LAT, + orbit, + ) + n_samples = 1e3 + arange_len = int(np.ceil(n_samples / (orbit.Np * orbit.Nsp))) + candidate_times = (possible_ts[None, :] + np.arange(0, arange_len)[:, None] * orbit.orbital_period_sec).flatten() + + pos_vec = orbit.get_orbit_positions_time_instant( + candidate_times + ) + all_lats = pos_vec["lat"] + all_lons = pos_vec["lon"] + lat_err = np.abs(all_lats - GROUND_STA_LAT) + lon_err = np.abs(all_lons - GROUND_STA_LON) + # do not let us choose a satellite with lat_err > 0.1 + lon_err[lat_err > 0.1] = 100 + # you may choose other parameters, such as choosing a specific sat + # my_sat_i = 2 * orbit.Nsp + 1 + # my_sat_err = lon_err[my_sat_i].copy() + # lon_err[::] = 100 + # lon_err[my_sat_i] = my_sat_err + flat_selected_i = np.argmin(lon_err) + + selected_sat_i = flat_selected_i // len(candidate_times) + selected_t_i = flat_selected_i % len(candidate_times) + + sat_lat = all_lats[selected_sat_i][selected_t_i] + sat_lon = all_lons[selected_sat_i][selected_t_i] + selected_t = candidate_times[selected_t_i] + + print(f"Selected satellite {selected_sat_i} at time {selected_t}") + print(f"\tPositioned at (lat, lon) = ({sat_lat}, {sat_lon})") + print( + "\tAt an error of (dlat, dlon) =", + f"({np.abs(sat_lat - GROUND_STA_LAT)}, {np.abs(sat_lon - GROUND_STA_LON)})" + ) + # range = 20 + # range = orbit.orbital_period_sec + # t_start, t_end = selected_t - range, selected_t + range + t_start, t_end = min_and_max_t_from( + selected_t, + selected_sat_i, + orbit, + GROUND_STA_LAT, + GROUND_STA_LON, + GROUND_STA_ALT, + MIN_ELEV_ANGLE_DEG, + ) + print(f"Can simulate with t in [{t_start}, {t_end}]") + + orbit.t_min = t_start + orbit.t_max = t_end + + plot_orbit_trace(orbit).show() + plot_sampling( + orbit, + GROUND_STA_LAT, + GROUND_STA_LON, + GROUND_STA_ALT, + MIN_ELEV_ANGLE_DEG, + num_drops=1 + ).show() + + print("Orbit additional params:") + for k in [ + "Mo", "model_time_as_random_variable", "t_min", "t_max", + ]: + print(4 * " " + f"{k}:", getattr(orbit, k)) diff --git a/sharc/satellite/utils/sat_utils.py b/sharc/satellite/utils/sat_utils.py index a005aa9a0..49d86ed23 100644 --- a/sharc/satellite/utils/sat_utils.py +++ b/sharc/satellite/utils/sat_utils.py @@ -115,6 +115,25 @@ def calc_elevation(Le: np.ndarray, return np.degrees(elev_angle) +def haversine( + lon1: np.ndarray, + lat1: np.ndarray, + lon2: np.ndarray, + lat2: np.ndarray, + R=EARTH_RADIUS_M +): + """Calculates great-circle distance between 2 points on the surface of Earth. + Considers spherical earth. + + Returns np.ndarray with N distances in meters + """ + lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2]) + dlon = lon2 - lon1 + dlat = lat2 - lat1 + a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2 + return 2 * R * np.arcsin(np.sqrt(a)) + + if __name__ == "__main__": r1 = ecef2lla(7792.1450, 0, 0) print(r1) diff --git a/sharc/simulation.py b/sharc/simulation.py index a89bd627f..120807227 100644 --- a/sharc/simulation.py +++ b/sharc/simulation.py @@ -20,6 +20,7 @@ from sharc.station_manager import StationManager from sharc.results import Results from sharc.propagation.propagation_factory import PropagationFactory +from sharc.support.sharc_utils import wrap2_180, clip_angle class Simulation(ABC, Observable): @@ -482,6 +483,8 @@ def select_ue(self, random_number_gen: np.random.RandomState): self.ue, ) bs_active = np.where(self.bs.active)[0] + + assert np.all((-180 <= self.bs.azimuth) & (self.bs.azimuth <= 180)), "BS azimuth angles should be in [-180, 180] range" for bs in bs_active: # select K UE's among the ones that are connected to BS random_number_gen.shuffle(self.link[bs]) @@ -494,9 +497,14 @@ def select_ue(self, random_number_gen: np.random.RandomState): # add beam to BS antennas # limit beamforming angle - bs_beam_phi = np.clip( + beam_h_min, beam_h_max = wrap2_180( + self.parameters.imt.bs.antenna.array.horizontal_beamsteering_range + self.bs.azimuth[bs] + ) + + bs_beam_phi = clip_angle( self.bs_to_ue_phi[bs, ue], - *(self.parameters.imt.bs.antenna.array.horizontal_beamsteering_range + self.bs.azimuth[bs]) + beam_h_min, + beam_h_max, ) bs_beam_theta = np.clip( diff --git a/sharc/simulation_downlink.py b/sharc/simulation_downlink.py index 8555ca511..bbe977a9f 100644 --- a/sharc/simulation_downlink.py +++ b/sharc/simulation_downlink.py @@ -233,11 +233,17 @@ def calculate_sinr_ext(self): if self.overlapping_bandwidth > 0: # in_band_interf_power = self.param_system.tx_power_density + \ # 10 * np.log10(self.overlapping_bandwidth * 1e6) + 30 - in_band_interf_power = \ - self.param_system.tx_power_density + 10 * np.log10( - self.ue.bandwidth[ue, np.newaxis] * 1e6 - ) + 10 * np.log10(weights)[:, np.newaxis] - \ - self.coupling_loss_imt_system[ue, :][:, active_sys] + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + category=RuntimeWarning, + message="divide by zero encountered in log10", + ) + in_band_interf_power = \ + self.param_system.tx_power_density + 10 * np.log10( + self.ue.bandwidth[ue, np.newaxis] * 1e6 + ) + 10 * np.log10(weights)[:, np.newaxis] - \ + self.coupling_loss_imt_system[ue, :][:, active_sys] oob_power = np.resize(-500., (len(ue), 1)) if self.adjacent_channel: @@ -414,10 +420,105 @@ def calculate_external_interference(self): # applying a bandwidth scaling factor since UE transmits on a portion # of the interfered systems bandwidth # calculate interference only from active UE's - rx_interference = 0 + pow_coch = -np.inf + # These are in dB. Turn to zero linear. + tx_oob = -np.inf + rx_oob = -np.inf bs_active = np.where(self.bs.active)[0] + # this implm assumes some parameters will be same for all interferring BS's + frst_bs = bs_active[0] + if self.co_channel: + ue = self.link[frst_bs] + weights = self.calculate_bw_weights( + self.ue.bandwidth[ue], + self.ue.center_freq[ue], + self.param_system.bandwidth, + 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": + # The unwanted emission is calculated in terms of TRP (after + # antenna). In SHARC implementation, ohmic losses are already + # included in coupling loss. Then, care has to be taken; + # otherwise ohmic loss will be included twice. + tx_oob = self.bs.spectral_mask.power_calc( + self.param_system.frequency, + self.system.bandwidth) + self.parameters.imt.bs.ohmic_loss + + elif self.parameters.imt.adjacent_ch_emissions == "ACLR": + non_overlap_sys_bw = self.param_system.bandwidth - self.overlapping_bandwidth + # NOTE: approximated equal to IMT bw + measurement_bw = self.parameters.imt.bandwidth + aclr_dB = self.parameters.imt.bs.adjacent_ch_leak_ratio + + if non_overlap_sys_bw > measurement_bw: + # NOTE: ACLR defines total leaked power over a fixed measurement bandwidth. + # If the victim bandwidth is wider, you’re assuming the same leakage + # profile extends beyond the ACLR-defined region, which may overestimate interference + # FIXME: if the victim bw fully contains tx bw, then + # EACH region should be <= measurement_bw + warn( + "Using IMT ACLR into system, but ACLR measurement bw is " + f"{measurement_bw} while the system bw is bigger ({non_overlap_sys_bw}).\n" + "Are you sure you intend to apply ACLR to the entire system bw?" + ) + + # tx_oob_in_measurement = (tx_pow_lin / aclr) + # => approx. PSD = (tx_pow_lin / aclr) / measurement_bw + # approximated received tx_oob = PSD * non_overlap_sys_bw + # NOTE: we don't get total power, but power per beam + # because later broadcast will sum this tx_oob `k` times + tx_oob = self.bs.tx_power[frst_bs] - aclr_dB + 10 * np.log10( + non_overlap_sys_bw / measurement_bw + ) + elif self.parameters.imt.adjacent_ch_emissions == "OFF": + pass + else: + raise ValueError( + f"No implementation for self.parameters.imt.adjacent_ch_emissions == { + self.parameters.imt.adjacent_ch_emissions}") + + # Calculate how much power is received in the adjacent channel + if self.param_system.adjacent_ch_reception == "ACS": + non_overlap_imt_bw = self.parameters.imt.bandwidth - self.overlapping_bandwidth + tx_bw = self.parameters.imt.bandwidth + acs_dB = self.param_system.adjacent_ch_selectivity + + # NOTE: only the power not overlapping is attenuated by ACS + # PSD = tx_pow_lin / tx_bw + # tx_pow_adj_lin = PSD * non_overlap_imt_bw + # rx_oob = tx_pow_adj_lin / acs + rx_oob = self.bs.tx_power[frst_bs] + 10 * np.log10( + non_overlap_imt_bw / tx_bw + ) - acs_dB + elif self.param_system.adjacent_ch_reception == "OFF": + if self.parameters.imt.adjacent_ch_emissions == "OFF": + raise ValueError( + "parameters.imt.adjacent_ch_emissions and parameters.imt.adjacent_ch_reception" + " cannot be both set to \"OFF\"") + else: + raise ValueError( + f"No implementation for self.param_system.adjacent_ch_reception == { + self.param_system.adjacent_ch_reception}") + sys_active = np.where(self.system.active)[0] + if len(sys_active) > 1: + raise NotImplementedError( + "Implementation does not support victim system with more than 1 active station" + ) + + rx_interference = 0 for bs in bs_active: active_beams = [ i for i in range( @@ -426,95 +527,12 @@ def calculate_external_interference(self): self.parameters.imt.ue.k, ) ] - if self.co_channel: - ue = self.link[bs] - weights = self.calculate_bw_weights( - self.ue.bandwidth[ue], - self.ue.center_freq[ue], - self.param_system.bandwidth, - self.param_system.frequency, - ) - - interference = self.bs.tx_power[bs] - \ - self.coupling_loss_imt_system[active_beams, sys_active] rx_interference += np.sum( - weights * np.power( - 10, - 0.1 * interference, - ), + 10 ** (0.1 * (pow_coch - self.coupling_loss_imt_system[active_beams, sys_active])) ) if self.adjacent_channel: - # These are in dB. Turn to zero linear. - tx_oob = -np.inf - rx_oob = -np.inf - # Calculate how much power is emitted in the adjacent channel: - if self.parameters.imt.adjacent_ch_emissions == "SPECTRAL_MASK": - # The unwanted emission is calculated in terms of TRP (after - # antenna). In SHARC implementation, ohmic losses are already - # included in coupling loss. Then, care has to be taken; - # otherwise ohmic loss will be included twice. - tx_oob = self.bs.spectral_mask.power_calc( - self.param_system.frequency, - self.system.bandwidth) + self.parameters.imt.bs.ohmic_loss - - elif self.parameters.imt.adjacent_ch_emissions == "ACLR": - non_overlap_sys_bw = self.param_system.bandwidth - self.overlapping_bandwidth - # NOTE: approximated equal to IMT bw - measurement_bw = self.parameters.imt.bandwidth - aclr_dB = self.parameters.imt.bs.adjacent_ch_leak_ratio - - if non_overlap_sys_bw > measurement_bw: - # NOTE: ACLR defines total leaked power over a fixed measurement bandwidth. - # If the victim bandwidth is wider, you’re assuming the same leakage - # profile extends beyond the ACLR-defined region, which may overestimate interference - # FIXME: if the victim bw fully contains tx bw, then - # EACH region should be <= measurement_bw - warn( - "Using IMT ACLR into system, but ACLR measurement bw is " - f"{measurement_bw} while the system bw is bigger ({non_overlap_sys_bw}).\n" - "Are you sure you intend to apply ACLR to the entire system bw?" - ) - - # tx_oob_in_measurement = (tx_pow_lin / aclr) - # => approx. PSD = (tx_pow_lin / aclr) / measurement_bw - # approximated received tx_oob = PSD * non_overlap_sys_bw - # NOTE: we don't get total power, but power per beam - # because later broadcast will sum this tx_oob `k` times - tx_oob = self.bs.tx_power[bs] - aclr_dB + 10 * np.log10( - non_overlap_sys_bw / measurement_bw - ) - elif self.parameters.imt.adjacent_ch_emissions == "OFF": - pass - else: - raise ValueError( - f"No implementation for self.parameters.imt.adjacent_ch_emissions == { - self.parameters.imt.adjacent_ch_emissions}") - - # Calculate how much power is received in the adjacent channel - if self.param_system.adjacent_ch_reception == "ACS": - non_overlap_imt_bw = self.parameters.imt.bandwidth - self.overlapping_bandwidth - tx_bw = self.parameters.imt.bandwidth - acs_dB = self.param_system.adjacent_ch_selectivity - - # NOTE: only the power not overlapping is attenuated by ACS - # PSD = tx_pow_lin / tx_bw - # tx_pow_adj_lin = PSD * non_overlap_imt_bw - # rx_oob = tx_pow_adj_lin / acs - rx_oob = self.bs.tx_power[bs] + 10 * np.log10( - non_overlap_imt_bw / tx_bw - ) - acs_dB - elif self.param_system.adjacent_ch_reception == "OFF": - if self.parameters.imt.adjacent_ch_emissions == "OFF": - raise ValueError( - "parameters.imt.adjacent_ch_emissions and parameters.imt.adjacent_ch_reception" - " cannot be both set to \"OFF\"") - else: - raise ValueError( - f"No implementation for self.param_system.adjacent_ch_reception == { - self.param_system.adjacent_ch_reception}") - # oob_power per beam # NOTE: we only consider one beam since all beams should have gain # of a single element for IMT, and as such the coupling loss should be the @@ -526,15 +544,17 @@ def calculate_external_interference(self): # so more would have to be fixed before this assert np.all(adj_loss == adj_loss.flat[0]) - tx_oob -= adj_loss[0, :] + tx_oob_s = tx_oob - adj_loss[0, :] if self.param_system.adjacent_ch_reception != "OFF": - rx_oob -= self.coupling_loss_imt_system[active_beams, sys_active] + rx_oob_s = rx_oob - self.coupling_loss_imt_system[active_beams, sys_active] + else: + rx_oob_s = -np.inf # Out of band power # sum linearly power leaked into band and power received in the # adjacent band oob_power = 10 * np.log10( - 10 ** (0.1 * tx_oob) + 10 ** (0.1 * rx_oob) + 10 ** (0.1 * tx_oob_s) + 10 ** (0.1 * rx_oob_s) ) # System rx interference diff --git a/sharc/simulation_uplink.py b/sharc/simulation_uplink.py index 37c4db1c5..2a4e48bbb 100644 --- a/sharc/simulation_uplink.py +++ b/sharc/simulation_uplink.py @@ -558,9 +558,9 @@ def collect_results(self, write_to_file: bool, snapshot_number: int): self.results.imt_system_antenna_gain.extend( self.imt_system_antenna_gain[np.ix_(sys_active, active_beams)].flatten(), ) - self.results.imt_system_antenna_gain_adjacent.extend( - self.imt_system_antenna_gain_adjacent[np.ix_(sys_active, active_beams)].flatten(), - ) + if len(self.imt_system_antenna_gain_adjacent): + self.results.imt_system_antenna_gain_adjacent.extend( + self.imt_system_antenna_gain_adjacent[np.ix_(sys_active, active_beams)].flatten(),) self.results.imt_system_path_loss.extend( self.imt_system_path_loss[np.ix_(sys_active, active_beams)].flatten(), ) diff --git a/sharc/station_factory.py b/sharc/station_factory.py index f6c2dc0d8..22c0bae14 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 GeometryConverter +from sharc.support.sharc_utils import wrap2_180 class StationFactory(object): @@ -119,7 +120,7 @@ def generate_imt_base_stations( else: imt_base_stations.height = param.bs.height * np.ones(num_bs) - imt_base_stations.azimuth = topology.azimuth + imt_base_stations.azimuth = wrap2_180(topology.azimuth) imt_base_stations.active = random_number_gen.rand( num_bs, ) < param.bs.load_probability @@ -154,11 +155,12 @@ def generate_imt_base_stations( num_bs, dtype=Antenna, ) - for i in range(num_bs): - imt_base_stations.antenna[i] = \ - AntennaFactory.create_antenna( - param.bs.antenna, imt_base_stations.azimuth[i], - imt_base_stations.elevation[i],) + imt_base_stations.antenna = AntennaFactory.create_n_antennas( + param.bs.antenna, + imt_base_stations.azimuth, + imt_base_stations.elevation, + num_bs + ) # imt_base_stations.antenna = [AntennaOmni(0) for bs in range(num_bs)] imt_base_stations.bandwidth = param.bandwidth * np.ones(num_bs) @@ -457,11 +459,12 @@ def generate_imt_ue_outdoor( # TODO: this piece of code works only for uplink ue_param_ant.get_antenna_parameters() - for i in range(num_ue): - imt_ue.antenna[i] = AntennaFactory.create_antenna( - param.ue.antenna, imt_ue.azimuth[i], - imt_ue.elevation[i], - ) + imt_ue.antenna = AntennaFactory.create_n_antennas( + param.ue.antenna, + imt_ue.azimuth, + imt_ue.elevation, + num_ue, + ) # imt_ue.antenna = [AntennaOmni(0) for bs in range(num_ue)] imt_ue.bandwidth = param.bandwidth * np.ones(num_ue) @@ -639,7 +642,7 @@ def generate_imt_ue_indoor( for i in range(num_ue): imt_ue.antenna[i] = AntennaBeamformingImt( par, imt_ue.azimuth[i], - imt_ue.elevation[i], ue_param_ant.subarray + imt_ue.elevation[i], ) # imt_ue.antenna = [AntennaOmni(0) for bs in range(num_ue)] @@ -844,6 +847,7 @@ def generate_single_space_station( space_station.elevation[0]) ]) + space_station.z = space_station.height space_station.bandwidth = param.bandwidth space_station.noise_temperature = param.noise_temperature space_station.thermal_noise = -500 @@ -1188,37 +1192,11 @@ def generate_single_earth_station( [param.geometry.elevation.fixed], ) - match param.antenna.pattern: - case "OMNI": - single_earth_station.antenna = np.array( - [AntennaOmni(param.antenna.gain)], - ) - case "ITU-R S.465": - single_earth_station.antenna = np.array( - [AntennaS465(param.antenna.itu_r_s_465)], - ) - case "ITU-R Reg. RR. Appendice 7 Annex 3": - single_earth_station.antenna = np.array( - [AntennaReg_RR_A7_3(param.antenna.itu_reg_rr_a7_3)], - ) - case "ITU-R S.1855": - single_earth_station.antenna = np.array( - [AntennaS1855(param.antenna.itu_r_s_1855)], - ) - case "MODIFIED ITU-R S.465": - single_earth_station.antenna = np.array( - [AntennaModifiedS465(param.antenna.itu_r_s_465_modified)], - ) - case "ITU-R S.580": - single_earth_station.antenna = np.array( - [AntennaS580(param.antenna.itu_r_s_580)], - ) - case _: - sys.stderr.write( - "ERROR\nInvalid FSS ES antenna pattern: " + - param.antenna_pattern, - ) - sys.exit(1) + single_earth_station.antenna = np.array([ + AntennaFactory.create_antenna( + param.antenna, single_earth_station.azimuth, single_earth_station.elevation + ) + ]) single_earth_station.active = np.array([True]) single_earth_station.bandwidth = np.array([param.bandwidth]) @@ -1247,7 +1225,9 @@ def generate_single_earth_station( else: raise ValueError(f"Invalid or not implemented spectral mask - {param.spectral_mask}") - single_earth_station.spectral_mask.set_mask(param.tx_power_density + 10 * np.log10(param.bandwidth * 1e6)) + single_earth_station.spectral_mask.set_mask( + param.tx_power_density + 10 * np.log10(param.bandwidth * 1e6) + 30 + ) return single_earth_station @@ -1630,7 +1610,8 @@ def generate_mss_ss(param_mss: ParametersMssSs): 10 * np.log10( param_mss.bandwidth * - 1e6)) + 1e6) + 30 + ) return mss_ss @@ -1695,7 +1676,8 @@ def generate_mss_d2d( 10 * np.log10( params.bandwidth * - 1e6)) + 1e6) + 30 + ) # Configure satellite positions in the StationManager mss_d2d.x = mss_d2d_values["sat_x"] diff --git a/sharc/station_manager.py b/sharc/station_manager.py index 81dc6c382..9ab7a9088 100644 --- a/sharc/station_manager.py +++ b/sharc/station_manager.py @@ -149,14 +149,17 @@ def get_3d_distance_to(self, station) -> np.array: np.array 3D distance matrix between stations. """ - distance = np.empty([self.num_stations, station.num_stations]) - for i in range(self.num_stations): - distance[i] = np.sqrt( - np.power(self.x[i] - station.x, 2) + - np.power(self.y[i] - station.y, 2) + - np.power(self.z[i] - station.z, 2) - ) - return distance + dx = np.subtract.outer(self.x, station.x).astype(np.float64) + dy = np.subtract.outer(self.y, station.y).astype(np.float64) + dz = np.subtract.outer(self.z, station.z).astype(np.float64) + np.square(dx, out=dx) + np.square(dy, out=dy) + np.square(dz, out=dz) + np.sqrt( + dx + dy + dz, + out=dx + ) + return dx def get_dist_angles_wrap_around(self, station) -> np.array: """Calculate distances and angles using the wrap-around technique. @@ -286,20 +289,22 @@ def get_pointing_vector_to(self, station) -> tuple: theta is calculated with respect to z counter-clockwise). """ - point_vec_x = station.x - self.x[:, np.newaxis] - point_vec_y = station.y - self.y[:, np.newaxis] - point_vec_z = station.z - self.z[:, np.newaxis] + # malloc + dx = (station.x - self.x[:, np.newaxis]).astype(np.float64) + dy = (station.y - self.y[:, np.newaxis]).astype(np.float64) + dz = (station.z - self.z[:, np.newaxis]).astype(np.float64) dist = self.get_3d_distance_to(station) - phi = np.array( - np.rad2deg( - np.arctan2( - point_vec_y, point_vec_x, - ), - ), ndmin=2, - ) - theta = np.rad2deg(np.arccos(point_vec_z / dist)) + # NOTE: doing in place calculations + phi = np.rad2deg(np.arctan2(dy, dx, out=dx), out=dx) + # delete reference dx + del dx + + # in place calculations + theta = np.rad2deg(np.arccos(np.clip(dz / dist, -1.0, 1.0, out=dz), out=dz), out=dz) + # delete reference dz + del dz return phi, theta diff --git a/sharc/support/named_tuples.py b/sharc/support/named_tuples.py deleted file mode 100644 index 75945ad5c..000000000 --- a/sharc/support/named_tuples.py +++ /dev/null @@ -1,16 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Thu May 11 14:21:57 2017 - -@author: Calil -""" - -from collections import namedtuple - -AntennaPar = namedtuple( - "AntennaPar", - "adjacent_antenna_model normalization normalization_data element_pattern \ - element_max_g element_phi_3db element_theta_3db element_am element_sla_v n_rows \ - n_columns element_horiz_spacing element_vert_spacing multiplication_factor \ - minimum_array_gain downtilt", -) diff --git a/sharc/support/sharc_geom.py b/sharc/support/sharc_geom.py index da00032d8..54317dae4 100644 --- a/sharc/support/sharc_geom.py +++ b/sharc/support/sharc_geom.py @@ -608,6 +608,8 @@ def shrink_countries_by_km( def generate_grid_in_polygon( polygon: shp.geometry.Polygon, hexagon_radius: float, + rotation_deg: typing.Optional[float] = None, + translation: typing.Optional[typing.Tuple[float, float]] = None, ): """Generate a hexagonal grid inside a polygon and return points in EARTH_DEFAULT_CRS. @@ -621,6 +623,11 @@ def generate_grid_in_polygon( Polygon to fill with a grid. hexagon_radius : float Radius of the hexagons in the grid (in meters). + rotation_deg : float or None, optional + Rotation of the grid in degrees (default: None, no rotation). + translation : tuple of float or None, optional + Translation of the grid in meters (dx, dy) (default: None, no translation). + The translation must not exceed the grid spacing in magnitude. Returns ------- @@ -644,14 +651,18 @@ def generate_grid_in_polygon( # Transform to projection where unit is meters polygon_proj = shp.ops.transform(to_proj, polygon) - # create a bounding box, afterwards we filter to polygon site - minx, miny, maxx, maxy = polygon_proj.bounds + # Determine x/y spacing x_spacing = 3 * hexagon_radius y_spacing = hexagon_radius * np.sqrt(3) / 2 - x_vals = np.arange(minx, maxx + x_spacing, x_spacing) + # Create a bounding circle, afterwards we filter to polygon site + cx, cy = polygon_proj.centroid.coords[0] + minx, miny, maxx, maxy = polygon_proj.bounds + bbox_diag = np.hypot(maxx - minx, maxy - miny) + bound_radius = bbox_diag / 2 + 3 * hexagon_radius - y_vals = np.arange(miny, maxy + y_spacing, y_spacing) + x_vals = np.arange(cx - bound_radius, cx + bound_radius + x_spacing, x_spacing) + y_vals = np.arange(cy - bound_radius, cy + bound_radius + y_spacing, y_spacing) x_vals, y_vals = np.meshgrid(x_vals, y_vals) @@ -661,6 +672,20 @@ def generate_grid_in_polygon( x_vals = x_vals.ravel() y_vals = y_vals.ravel() + if rotation_deg is not None: + theta = np.deg2rad(rotation_deg) + cos_t, sin_t = np.cos(theta), np.sin(theta) + x_rot = cos_t * (x_vals - cx) - sin_t * (y_vals - cy) + cx + y_rot = sin_t * (x_vals - cx) + cos_t * (y_vals - cy) + cy + x_vals, y_vals = x_rot, y_rot + + if translation is not None: + dx, dy = translation + if abs(dx) > x_spacing or abs(dy) > y_spacing: + raise ValueError("generate_grid_in_polygon.translation (dx, dy) must not exceed grid spacing (x_spacing, y_spacing) in magnitude") + x_vals += dx + y_vals += dy + # Return to EARTH_DEFAULT_CRS xt, yt = from_proj(x_vals, y_vals) @@ -676,34 +701,66 @@ def generate_grid_in_polygon( def generate_grid_in_multipolygon( poly: typing.Union[shp.geometry.MultiPolygon, shp.geometry.Polygon], - km: float + km: float, + random_transform_on_grid: bool = False, + rng: np.random.RandomState = None, ) -> list[shp.geometry.MultiPolygon]: - """Generate a grid in a MultiPolygon or Polygon, shrinking each by a given number of kilometers. + """Generate a hexagonal grid in a MultiPolygon or Polygon, + considering a hexagon radius in km. - For each polygon, create a grid and return a single 2xN array of longitudes and latitudes. + For each polygon, create a grid and return a single 2xN array of longitudes and latitudes + containing all grids. Parameters ---------- poly : typing.Union[shp.geometry.MultiPolygon, shp.geometry.Polygon] The MultiPolygon or Polygon to process. km : float - Number of kilometers to shrink each polygon by. + Hexagon radius in km + random_transform_on_grid : bool, optional + Whether to apply a random rotation and translation to the grid (default: False). + rng : np.random.RandomState, optional + Random number generator to use if random_transform_on_grid is True (default: None). Returns ------- np.ndarray 2xN array: first row is longitudes, second row is latitudes. """ + if random_transform_on_grid: + assert rng is not None + lons = [] lats = [] if poly.geom_type == 'Polygon': - x, y = generate_grid_in_polygon(poly, km) + if random_transform_on_grid: + x, y = generate_grid_in_polygon( + poly, km, + rng.uniform(-180., 180.), + ( + 3 * rng.uniform(-km, km), + rng.uniform(-km, km) * np.sqrt(3) / 2 + ) + ) + else: + x, y = generate_grid_in_polygon(poly, km) lons.extend(x) lats.extend(y) elif poly.geom_type == 'MultiPolygon': for p in poly.geoms: - x, y = generate_grid_in_polygon(p, km) + if random_transform_on_grid: + x, y = generate_grid_in_polygon( + p, km, + rng.uniform(-180., 180.), + ( + 3 * rng.uniform(-km, km), + rng.uniform(-km, km) * np.sqrt(3) / 2 + ) + ) + else: + x, y = generate_grid_in_polygon(p, km) + lons.extend(x) lats.extend(y) diff --git a/sharc/support/sharc_utils.py b/sharc/support/sharc_utils.py index bc7644d53..c5d12398d 100644 --- a/sharc/support/sharc_utils.py +++ b/sharc/support/sharc_utils.py @@ -26,6 +26,49 @@ def is_float(s: str) -> bool: return False +def wrap2_180(a): + """ + Wraps angles to the [-180, 180) range + """ + return (a + 180) % 360 - 180 + + +def angular_dist(a1, a2): + """ + Returns smallest angular distance between 2 angles + in the [0, 180] range + """ + return np.abs(wrap2_180(a1 - a2)) + + +def clip_angle(a, a_min, a_max): + """ + Returns `a` if it is inside the angle range defined by [`a_min`, `a_max`] + otherwise returns the closest bound, be it `a_min` or `a_max` + It is assumed all angles passes are already in [-180, 180] range + """ + outside_rng = False + # NOTE: angle limits such as -180 + [-60, 60] + # should be passed as wrapped around [120, -120] + # So it is important to check for this case + if a_min > a_max: + if a > a_max and a < a_min: + outside_rng = True + else: + # Normal interval + if a < a_min or a > a_max: + outside_rng = True + + # always clip to the closest bound + if outside_rng: + if angular_dist(a, a_min) < angular_dist(a, a_max): + return a_min + else: + return a_max + + return a + + def to_scalar(x): """Convert a numpy scalar or array to a Python scalar if possible.""" if isinstance(x, np.ndarray): diff --git a/sharc/topology/topology_factory.py b/sharc/topology/topology_factory.py index 75bad88e5..646df3934 100644 --- a/sharc/topology/topology_factory.py +++ b/sharc/topology/topology_factory.py @@ -27,7 +27,8 @@ def createTopology(parameters: Parameters, if parameters.imt.topology.type == "SINGLE_BS": return TopologySingleBaseStation( parameters.imt.topology.single_bs.cell_radius, - parameters.imt.topology.single_bs.num_clusters + parameters.imt.topology.single_bs.num_clusters, + parameters.imt.topology.single_bs.azimuth ) elif parameters.imt.topology.type == "MACROCELL": return TopologyMacrocell( diff --git a/sharc/topology/topology_imt_mss_dc.py b/sharc/topology/topology_imt_mss_dc.py index f8f14266e..292f37b5f 100644 --- a/sharc/topology/topology_imt_mss_dc.py +++ b/sharc/topology/topology_imt_mss_dc.py @@ -65,23 +65,6 @@ def __init__(self, params: ParametersImtMssDc, self.lat = None self.lon = None - self.orbits = [] - - # Iterate through each orbit defined in the parameters - for param in self.orbit_params.orbits: - # Instantiate an OrbitModel for the current orbit - 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 - ) - self.orbits.append(orbit) @staticmethod def get_coordinates( @@ -139,7 +122,11 @@ def get_coordinates( 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 + 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( @@ -386,7 +373,10 @@ def get_satellite_pointing( (before sharc's coordinate transformation, with distances in km) """ if orbit_params.beam_positioning.type == "SERVICE_GRID": - orbit_params.beam_positioning.service_grid.validate("service_grid") + # TODO: remove this reset_grid call from here. + # Since it should be called once per drop, it shouldn't + # be this deep in the program + orbit_params.beam_positioning.service_grid.reset_grid("service_grid", random_number_gen) # 2xN, (lon, lat) grid = orbit_params.beam_positioning.service_grid.lon_lat_grid diff --git a/sharc/topology/topology_single_base_station.py b/sharc/topology/topology_single_base_station.py index 13bc55273..a1719815d 100644 --- a/sharc/topology/topology_single_base_station.py +++ b/sharc/topology/topology_single_base_station.py @@ -5,6 +5,7 @@ @author: edgar """ import numpy as np +from typing import List import matplotlib.pyplot as plt import matplotlib.axes import matplotlib.patches as patches @@ -22,13 +23,18 @@ class TopologySingleBaseStation(Topology): AZIMUTH = [0, 180] ALLOWED_NUM_CLUSTERS = [1, 2] - def __init__(self, cell_radius: float, num_clusters: int): + def __init__(self, cell_radius: float, num_clusters: int, azimuth: List[float] | str | None = None): """ Constructor method that sets the object attributes Parameters ---------- cell_radius : radius of the cell + num_clusters : number of clusters (1 or 2) + azimuth : List of azimuths for each cluster in degrees or "random" for random azimuths + Raises + ------ + ValueError : if num_clusters is not 1 or 2 """ if num_clusters not in TopologySingleBaseStation.ALLOWED_NUM_CLUSTERS: error_message = "invalid number of clusters ({})".format( @@ -36,6 +42,10 @@ def __init__(self, cell_radius: float, num_clusters: int): ) raise ValueError(error_message) + if isinstance(azimuth, str): + if azimuth.lower() != "random": + raise ValueError("azimuth must be a list of floats or 'random'") + self.azimuth_param = azimuth intersite_distance = 2 * cell_radius super().__init__(intersite_distance, cell_radius) self.num_clusters = num_clusters @@ -49,16 +59,23 @@ def calculate_coordinates(self, random_number_gen=np.random.RandomState()): if self.num_clusters == 1: self.x = np.array([0]) self.y = np.array([0]) - self.azimuth = TopologySingleBaseStation.AZIMUTH[0] * np.ones( - 1, - ) self.num_base_stations = 1 elif self.num_clusters == 2: self.x = np.array([0, self.intersite_distance]) self.y = np.array([0, 0]) - self.azimuth = np.array(TopologySingleBaseStation.AZIMUTH) self.num_base_stations = 2 self.indoor = np.zeros(self.num_base_stations, dtype=bool) + # Set the azimuth accordingly + if self.azimuth_param is None: + self.azimuth = np.array(self.AZIMUTH[:self.num_clusters]).astype(float) + elif isinstance(self.azimuth_param, list): + if len(self.azimuth_param) != self.num_clusters: + raise ValueError( + "Length of azimuth list must be equal to num_clusters") + self.azimuth = np.array(self.azimuth_param).astype(float) + elif self.azimuth_param.lower() == "random": + self.azimuth = random_number_gen.randint(0, 360, self.num_clusters).astype(float) + self.z = np.zeros_like(self.x) def plot(self, ax: matplotlib.axes.Axes): @@ -87,7 +104,7 @@ def plot(self, ax: matplotlib.axes.Axes): if __name__ == '__main__': cell_radius = 100 - num_clusters = 2 + num_clusters = 1 topology = TopologySingleBaseStation(cell_radius, num_clusters) topology.calculate_coordinates() diff --git a/tests/_test_beamforming_normalizer.py b/tests/_test_beamforming_normalizer.py index e607b1f7e..5adba1d9e 100644 --- a/tests/_test_beamforming_normalizer.py +++ b/tests/_test_beamforming_normalizer.py @@ -11,7 +11,7 @@ import os from sharc.antenna.beamforming_normalization.beamforming_normalizer import BeamformingNormalizer -from sharc.support.named_tuples import AntennaPar +from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt from sharc.antenna.antenna_beamforming_imt import AntennaBeamformingImt @@ -40,23 +40,25 @@ def setUp(self): multiplication_factor = 12 minimum_array_gain = -200 downtilt = 0 - self.par_1 = AntennaPar( - adjacent_antenna_model, - norm, - norm_file, - element_pattern, - element_max_g, - element_phi_3db, - element_theta_3db, - element_am, - element_sla_v, - n_rows, - n_columns, - horiz_spacing, - vert_spacing, - multiplication_factor, - minimum_array_gain, - downtilt, + self.par_1 = ParametersAntennaImt( + horizontal_beamsteering_range=[-180., 180.], + vertical_beamsteering_range=[-90., 90.], + adjacent_antenna_model=adjacent_antenna_model, + normalization=norm, + normalization_file=norm_file, + element_pattern=element_pattern, + element_max_g=element_max_g, + element_phi_3db=element_phi_3db, + element_theta_3db=element_theta_3db, + element_am=element_am, + element_sla_v=element_sla_v, + n_rows=n_rows, + n_columns=n_columns, + element_horiz_spacing=horiz_spacing, + element_vert_spacing=vert_spacing, + multiplication_factor=multiplication_factor, + minimum_array_gain=minimum_array_gain, + downtilt=downtilt, ) # Test 2: UE configuration @@ -68,8 +70,8 @@ def setUp(self): norm_file = None element_pattern = "M2101" element_max_g = 5 - element_phi_deg_3db = 90 - element_theta_deg_3db = 90 + element_phi_3_db = 90 + element_theta_3_db = 90 element_am = 25 element_sla_v = 25 n_rows = 4 @@ -79,23 +81,25 @@ def setUp(self): multiplication_factor = 12 minimum_array_gain = -200 downtilt = 0 - self.par_2 = AntennaPar( - adjacent_antenna_model, - norm, - norm_file, - element_pattern, - element_max_g, - element_phi_deg_3db, - element_theta_deg_3db, - element_am, - element_sla_v, - n_rows, - n_columns, - horiz_spacing, - vert_spacing, - multiplication_factor, - minimum_array_gain, - downtilt, + self.par_2 = ParametersAntennaImt( + horizontal_beamsteering_range=[-180., 180.], + vertical_beamsteering_range=[-90., 90.], + adjacent_antenna_model=adjacent_antenna_model, + normalization=norm, + normalization_file=norm_file, + element_pattern=element_pattern, + element_max_g=element_max_g, + element_phi_3db=element_phi_3_db, + element_theta_3db=element_theta_3_db, + element_am=element_am, + element_sla_v=element_sla_v, + n_rows=n_rows, + n_columns=n_columns, + element_horiz_spacing=horiz_spacing, + element_vert_spacing=vert_spacing, + multiplication_factor=multiplication_factor, + minimum_array_gain=minimum_array_gain, + downtilt=downtilt, ) # Test 3: BS configuration @@ -107,8 +111,8 @@ def setUp(self): norm_file = None element_pattern = "M2101" element_max_g = 5 - element_phi_deg_3db = 65 - element_theta_deg_3db = 65 + element_phi_3_db = 65 + element_theta_3_db = 65 element_am = 30 element_sla_v = 30 n_rows = 8 @@ -118,23 +122,25 @@ def setUp(self): multiplication_factor = 12 minimum_array_gain = -200 downtilt = 0 - self.par_3 = AntennaPar( - adjacent_antenna_model, - norm, - norm_file, - element_pattern, - element_max_g, - element_phi_deg_3db, - element_theta_deg_3db, - element_am, - element_sla_v, - n_rows, - n_columns, - horiz_spacing, - vert_spacing, - multiplication_factor, - minimum_array_gain, - downtilt, + self.par_3 = ParametersAntennaImt( + horizontal_beamsteering_range=[-180., 180.], + vertical_beamsteering_range=[-90., 90.], + adjacent_antenna_model=adjacent_antenna_model, + normalization=norm, + normalization_file=norm_file, + element_pattern=element_pattern, + element_max_g=element_max_g, + element_phi_3db=element_phi_3_db, + element_theta_3db=element_theta_3_db, + element_am=element_am, + element_sla_v=element_sla_v, + n_rows=n_rows, + n_columns=n_columns, + element_horiz_spacing=horiz_spacing, + element_vert_spacing=vert_spacing, + multiplication_factor=multiplication_factor, + minimum_array_gain=minimum_array_gain, + downtilt=downtilt, ) def test_construction(self): diff --git a/tests/parameters/parameters_for_testing.yaml b/tests/parameters/parameters_for_testing.yaml index a815020a5..20297abad 100644 --- a/tests/parameters/parameters_for_testing.yaml +++ b/tests/parameters/parameters_for_testing.yaml @@ -97,6 +97,15 @@ imt: # Number of clusters in single base station topology # You can simulate 1 or 2 BS's with this topology num_clusters: 2 + ########################################################################### + # Azimuth of each cluster in single base station topology [deg] + # You can provide a list of azimuths with length equal to num_clusters + # or use "random" to have them randomly generated + azimuth: + - 60.0 + - 240.0 + # azimuth: !!python/tuple [0.0, 180.0] + # azimuth: random ########################################################################### # Hotspot Topology parameters. Relevant when imt.topology.type == "HOTSPOT" hotspot: @@ -223,6 +232,7 @@ imt: min: 0.0 max: 66.1 service_grid: + transform_grid_randomly: true # add per drop rand rotation + transl. to grid # by default this already gets shapefile from natural earth country_shapes_filename: sharc/data/countries/ne_110m_admin_0_countries.shp # By default, these are the same as taken from `sat_is_active_if.lat_long_inside_country` @@ -279,6 +289,7 @@ imt: sats_per_plane: 20 long_asc_deg: 30.0 phasing_deg: 2.0 + model_time_as_random_variable: true - n_planes: 26 inclination_deg: 97.77 perigee_alt_km: 595 @@ -286,6 +297,9 @@ imt: sats_per_plane: 30 long_asc_deg: 14.0 phasing_deg: 7.8 + model_time_as_random_variable: true + t_min: 1.0 + t_max: 100.0 ########################################################################### # Defines the antenna model to be used in compatibility studies between # IMT and other services in adjacent band @@ -566,9 +580,6 @@ imt: shadowing : FALSE noise_temperature : 290 param_p619: - ########################################################################### - # altitude of the space station [m] - space_station_alt_m: 540000 ########################################################################### # altitude of ES system [m] earth_station_alt_m: 1200 @@ -576,9 +587,9 @@ imt: # latitude of ES system [m] earth_station_lat_deg: 13 ########################################################################### - # difference between longitudes of IMT-NTN station and FSS-ES system [deg] - # (positive if space-station is to the East of earth-station) - earth_station_long_diff_deg: 10 + mean_clutter_height: 'low' + ########################################################################### + below_rooftop: 50 haps: ########################################################################### # HAPS center frequency [MHz] @@ -752,9 +763,6 @@ single_earth_station: ########################################################################### # Parameters for P619 channel model param_p619: - ########################################################################### - # altitude of the space station [m] - space_station_alt_m: 540000 ########################################################################### # altitude of ES system [m] earth_station_alt_m: 1200 @@ -762,10 +770,9 @@ single_earth_station: # latitude of ES system [m] earth_station_lat_deg: 13 ########################################################################### - # difference between longitudes of IMT-NTN station and FSS-ES system [deg] - # (positive if space-station is to the East of earth-station) - earth_station_long_diff_deg: 10 - + mean_clutter_height: 'low' + ########################################################################### + below_rooftop: 50 ########################################################################### # season - season of the year. # Only actually useful for P619, but we can't put it inside param_p619 without changing more code @@ -813,6 +820,7 @@ single_earth_station: ########################################################################### # determine whether clutter loss following ITU-R P.2108 is added (TRUE/FALSE) clutter_loss: TRUE + clutter_type: 'one-end' # HDFSS propagation parameters param_hdfss: @@ -1001,6 +1009,7 @@ mss_d2d: min: 0.0 max: 66.1 service_grid: + transform_grid_randomly: true # add per drop rand rotation + transl. to grid # by default this already gets shapefile from natural earth country_shapes_filename: sharc/data/countries/ne_110m_admin_0_countries.shp # By default, these are the same as taken from `sat_is_active_if.lat_long_inside_country` @@ -1071,12 +1080,10 @@ mss_d2d: # latitude of ES system [m] earth_station_lat_deg: 0 ########################################################################### - # difference between longitudes of IMT-NTN station and FSS-ES system [deg] - # (positive if space-station is to the East of earth-station) - earth_station_long_diff_deg: 0.75 - ########################################################################### # year season: SUMMER of WINTER season: SUMMER + mean_clutter_height: 'low' + below_rooftop: 50 # Parameters for the orbits orbits: - n_planes: 20 @@ -1093,10 +1100,14 @@ mss_d2d: sats_per_plane: 20 long_asc_deg: 30.0 phasing_deg: 2.0 + model_time_as_random_variable: true - n_planes: 26 inclination_deg: 97.77 perigee_alt_km: 595 apogee_alt_km: 595.0 sats_per_plane: 30 long_asc_deg: 14.0 - phasing_deg: 7.8 \ No newline at end of file + phasing_deg: 7.8 + model_time_as_random_variable: true + t_min: 1.0 + t_max: 100.0 diff --git a/tests/parameters/test_parameters.py b/tests/parameters/test_parameters.py index f41ea6d73..697693d3f 100644 --- a/tests/parameters/test_parameters.py +++ b/tests/parameters/test_parameters.py @@ -63,9 +63,9 @@ def test_parameters_imt(self): self.assertEqual( self.parameters.imt.ue.antenna.array.horizontal_beamsteering_range, (-180., - 180.)) + 179.9999)) self.assertEqual( - self.parameters.imt.ue.antenna.array.vertical_beamsteering_range, (0., 180.)) + self.parameters.imt.ue.antenna.array.vertical_beamsteering_range, (0., 179.9999)) self.assertEqual(self.parameters.imt.ue.k, 3) self.assertEqual(self.parameters.imt.ue.k_m, 1) self.assertEqual(self.parameters.imt.ue.indoor_percent, 5.0) @@ -224,6 +224,9 @@ def test_parameters_imt(self): self.assertEqual( self.parameters.imt.topology.single_bs.num_clusters, 2) + self.assertEqual( + self.parameters.imt.topology.single_bs.azimuth, [60.0, 240.0]) + """Test ParametersIndoor """ self.assertEqual( @@ -277,6 +280,9 @@ def test_parameters_imt(self): self.assertEqual( self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.beam_radius, 19000) + self.assertEqual( + self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.transform_grid_randomly, + True) self.assertEqual( self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.grid_margin_from_border, 0.11) @@ -327,6 +333,7 @@ def test_parameters_imt(self): 'sats_per_plane': 32, 'long_asc_deg': 18.0, 'phasing_deg': 3.9, + "model_time_as_random_variable": False, }, { 'n_planes': 12, @@ -336,6 +343,10 @@ def test_parameters_imt(self): 'sats_per_plane': 20, 'long_asc_deg': 30.0, 'phasing_deg': 2.0, + "model_time_as_random_variable": True, + # checking defaults + "t_min": 0.0, + "t_max": None, }, { 'n_planes': 26, @@ -345,31 +356,16 @@ def test_parameters_imt(self): 'sats_per_plane': 30, 'long_asc_deg': 14.0, 'phasing_deg': 7.8, + "model_time_as_random_variable": True, + "t_min": 1.0, + "t_max": 100.0, }, ] - for i, orbit_params in enumerate( - self.parameters.imt.topology.mss_dc.orbits): - self.assertEqual( - orbit_params.n_planes, - expected_orbit_params[i]['n_planes']) - self.assertEqual( - orbit_params.inclination_deg, - expected_orbit_params[i]['inclination_deg']) - self.assertEqual( - orbit_params.perigee_alt_km, - expected_orbit_params[i]['perigee_alt_km']) - self.assertEqual( - orbit_params.apogee_alt_km, - expected_orbit_params[i]['apogee_alt_km']) - self.assertEqual( - orbit_params.sats_per_plane, - expected_orbit_params[i]['sats_per_plane']) - self.assertEqual( - orbit_params.long_asc_deg, - expected_orbit_params[i]['long_asc_deg']) - self.assertEqual( - orbit_params.phasing_deg, - expected_orbit_params[i]['phasing_deg']) + for i, orbit_params in enumerate(self.parameters.imt.topology.mss_dc.orbits): + for k in expected_orbit_params[i].keys(): + self.assertEqual( + getattr(orbit_params, k), + expected_orbit_params[i][k]) def test_imt_validation(self): """ @@ -512,16 +508,8 @@ def test_parameters_single_earth_station(self): self.parameters.single_earth_station.param_p619.earth_station_alt_m, 1200, ) - self.assertEqual( - self.parameters.single_earth_station.param_p619.space_station_alt_m, - 540000, - ) self.assertEqual( self.parameters.single_earth_station.param_p619.earth_station_lat_deg, 13, ) - self.assertEqual( - self.parameters.single_earth_station.param_p619.earth_station_long_diff_deg, - 10, - ) self.assertEqual( self.parameters.single_earth_station.param_p452.atmospheric_pressure, 1, ) @@ -641,12 +629,13 @@ def test_parametes_mss_d2d(self): self.parameters.mss_d2d.param_p619.earth_station_alt_m, 0.0) self.assertEqual( self.parameters.mss_d2d.param_p619.earth_station_lat_deg, 0.0) - self.assertEqual( - self.parameters.mss_d2d.param_p619.earth_station_long_diff_deg, 0.0) self.assertEqual( self.parameters.mss_d2d.beam_positioning.service_grid.beam_radius, 19001) + self.assertEqual( + self.parameters.mss_d2d.beam_positioning.service_grid.transform_grid_randomly, + True) self.assertEqual( self.parameters.mss_d2d.beam_positioning.service_grid.grid_margin_from_border, 0.11) @@ -694,6 +683,7 @@ def test_parametes_mss_d2d(self): 'sats_per_plane': 32, 'long_asc_deg': 18.0, 'phasing_deg': 3.9, + "model_time_as_random_variable": False, }, { 'n_planes': 12, @@ -703,6 +693,10 @@ def test_parametes_mss_d2d(self): 'sats_per_plane': 20, 'long_asc_deg': 30.0, 'phasing_deg': 2.0, + "model_time_as_random_variable": True, + # check defaults + "t_min": 0.0, + "t_max": None, }, { 'n_planes': 26, @@ -712,30 +706,16 @@ def test_parametes_mss_d2d(self): 'sats_per_plane': 30, 'long_asc_deg': 14.0, 'phasing_deg': 7.8, + "model_time_as_random_variable": True, + "t_min": 1.0, + "t_max": 100.0, }, ] for i, orbit_params in enumerate(self.parameters.mss_d2d.orbits): - self.assertEqual( - orbit_params.n_planes, - expected_orbit_params[i]['n_planes']) - self.assertEqual( - orbit_params.inclination_deg, - expected_orbit_params[i]['inclination_deg']) - self.assertEqual( - orbit_params.perigee_alt_km, - expected_orbit_params[i]['perigee_alt_km']) - self.assertEqual( - orbit_params.apogee_alt_km, - expected_orbit_params[i]['apogee_alt_km']) - self.assertEqual( - orbit_params.sats_per_plane, - expected_orbit_params[i]['sats_per_plane']) - self.assertEqual( - orbit_params.long_asc_deg, - expected_orbit_params[i]['long_asc_deg']) - self.assertEqual( - orbit_params.phasing_deg, - expected_orbit_params[i]['phasing_deg']) + for k in expected_orbit_params[i].keys(): + self.assertEqual( + getattr(orbit_params, k), + expected_orbit_params[i][k]) def test_parameters_single_space_station(self): """Test ParametersSinglespaceStation @@ -818,10 +798,6 @@ def test_parameters_single_space_station(self): self.parameters.single_space_station.param_p619.earth_station_alt_m, self.parameters.single_space_station.geometry.es_altitude, ) - self.assertEqual( - self.parameters.single_space_station.param_p619.space_station_alt_m, - self.parameters.single_space_station.geometry.altitude, - ) self.assertEqual( self.parameters.single_space_station.param_p619.earth_station_lat_deg, self.parameters.single_space_station.geometry.es_lat_deg, diff --git a/tests/test_antenna_beamforming_imt.py b/tests/test_antenna_beamforming_imt.py index 95db72612..0137ef816 100644 --- a/tests/test_antenna_beamforming_imt.py +++ b/tests/test_antenna_beamforming_imt.py @@ -11,7 +11,6 @@ from sharc.antenna.antenna_beamforming_imt import AntennaBeamformingImt from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt -from sharc.support.named_tuples import AntennaPar class AntennaBeamformingImtTest(unittest.TestCase): @@ -479,28 +478,28 @@ def test_normalization(self): minimum_array_gain = -200 multiplication_factor = 12 downtilt = 0 - par = AntennaPar( - adjacent_antenna_model, - normalization, - norm_data, - element_pattern, - element_max_g, - element_phi_3db, - element_theta_3db, - element_am, - element_sla_v, - n_rows, - n_columns, - horiz_spacing, - vert_spacing, - multiplication_factor, - minimum_array_gain, - downtilt, + par = ParametersAntennaImt( + adjacent_antenna_model=adjacent_antenna_model, + normalization=normalization, + element_pattern=element_pattern, + element_max_g=element_max_g, + element_phi_3db=element_phi_3db, + element_theta_3db=element_theta_3db, + element_am=element_am, + element_sla_v=element_sla_v, + n_rows=n_rows, + n_columns=n_columns, + element_horiz_spacing=horiz_spacing, + element_vert_spacing=vert_spacing, + multiplication_factor=multiplication_factor, + minimum_array_gain=minimum_array_gain, + downtilt=downtilt, ) + par.normalization_data = norm_data # Create antenna objects self.antenna3 = AntennaBeamformingImt(par, 0.0, 0.0) # Normalized - par = par._replace(normalization=False) + par.normalization = False self.antenna4 = AntennaBeamformingImt(par, 0.0, 0.0) # Unormalized # Test co-channel gain: no beam diff --git a/tests/test_orbit_model.py b/tests/test_orbit_model.py index 29b19739c..795a807ae 100644 --- a/tests/test_orbit_model.py +++ b/tests/test_orbit_model.py @@ -1,6 +1,7 @@ import unittest import numpy as np from sharc.satellite.ngso.orbit_model import OrbitModel +from sharc.support.sharc_utils import angular_dist class TestOrbitModel(unittest.TestCase): @@ -18,6 +19,9 @@ def setUp(self): hp=1414, ha=1414, Mo=0, + model_time_as_random_variable=False, + t_min=0.0, + t_max=None ) def test_initialization(self): @@ -41,7 +45,7 @@ def test_orbit_parameters(self): 6845.3519, places=4) self.assertAlmostEqual(self.orbit.sat_sep_angle_deg, 60.0) - self.assertAlmostEqual(self.orbit.orbital_plane_inclination, 45.0) + self.assertAlmostEqual(self.orbit.orbital_plane_spacing, 45.0) def test_mean_anomalies(self): """Test mean anomaly calculations and satellite phasing logic.""" @@ -84,6 +88,58 @@ def test_mean_anomalies(self): np.testing.assert_array_almost_equal( np.diff(ma_deg, axis=1), r, decimal=4) + def test_limiting_time_range_for_random_drops(self): + """ + Testing if using model_time_as_random_variable + and t_min, t_max works correctly + """ + # if we limit time max to period / frac + # and considering a circular orbit + # then only 2 * pi / frac is traveled at max per each satellit + frac = 36 + self.orbit.t_max = self.orbit.orbital_period_sec / frac + max_traveled_angular_dist = 2 * np.pi / frac + + self.orbit.model_time_as_random_variable = True + # NOTE: choosing a larger frac lets you choose less samples + n_samples = int(2e2) + positions = self.orbit.get_orbit_positions_random( + rng=np.random.RandomState(123), + n_samples=n_samples, + ) + all_lats = positions["lat"] + all_lons = positions["lon"] + all_alts = positions["alt"] + all_xs = positions["sx"] + all_ys = positions["sy"] + all_zs = positions["sz"] + + self.assertEqual(all_lons.shape, all_lats.shape) + self.assertEqual(all_lons.shape, all_xs.shape) + self.assertEqual(all_lons.shape, all_ys.shape) + self.assertEqual(all_lons.shape, all_zs.shape) + self.assertEqual(all_lons.shape, all_alts.shape) + self.assertEqual(all_lons.shape, (self.orbit.Np * self.orbit.Nsp, n_samples)) + + for sat_lats, sat_lons in zip(all_lats, all_lons): + d_lon = np.deg2rad(angular_dist( + sat_lons[:, np.newaxis], sat_lons[np.newaxis, :] + )) + + # exhaustively get angular distance between all samples + a = np.pi / 2 - np.deg2rad(sat_lats[:, np.newaxis]) + b = np.pi / 2 - np.deg2rad(sat_lats[np.newaxis, :]) + cos_phi = np.cos(a) * np.cos(b) \ + + np.sin(a) * np.sin(b) * np.cos(d_lon) + phi = np.max(np.arccos( + # imprecision may accumulate enough for numbers to be slightly out + # of arccos range + np.clip(cos_phi, -1., 1.) + )) + + self.assertLessEqual(phi, max_traveled_angular_dist) + self.assertGreaterEqual(phi, max_traveled_angular_dist * 0.9) + if __name__ == '__main__': unittest.main() diff --git a/tests/test_propagation_clutter.py b/tests/test_propagation_clutter.py index a772d8ce5..8ddf7cf71 100644 --- a/tests/test_propagation_clutter.py +++ b/tests/test_propagation_clutter.py @@ -13,24 +13,29 @@ def setUp(self): def test_spatial_clutter_loss(self): """Test spatial clutter loss for different elevations and location percentages.""" - frequency = np.array([27000]) # MHz - elevation = np.array([0, 45, 90]) - loc_percentage = np.array([0.1, 0.5, 0.9]) - distance = np.array([1000]) # meters, dummy value - + frequency = np.array([27000, 27000, 27000]) # MHz + elevation = np.array([10, 20, 30]) + loc_percentage = np.array([.1, .5, .9]) + distance = np.array([1000, 1000, 1000]) # meters, dummy value + earth_station_height = np.array([10, 10, 10]) + mean_clutter_height = 'high' + below_rooftop = 100 loss = self.clutter_loss.get_loss( distance=distance, frequency=frequency, elevation=elevation, loc_percentage=loc_percentage, station_type=StationType.FSS_SS, + earth_station_height=earth_station_height, + mean_clutter_height=mean_clutter_height, + below_rooftop=below_rooftop, ) # Check the shape of the output self.assertEqual(loss.shape, (3,)) # Check if loss decreases with increasing elevation - self.assertTrue(loss[0] >= loss[1] >= loss[2]) + self.assertTrue(loss[0] <= loss[1] <= loss[2]) def test_terrestrial_clutter_loss(self): """Test terrestrial clutter loss for different frequencies and distances.""" @@ -38,12 +43,13 @@ def test_terrestrial_clutter_loss(self): distance = np.array([500, 2000]) # meters # Using a single value for location percentage loc_percentage = np.array([0.5]) - + clutter_type = 'one_end' loss = self.clutter_loss.get_loss( frequency=frequency, distance=distance, loc_percentage=loc_percentage, station_type=StationType.IMT_BS, + clutter_type=clutter_type ) self.assertEqual(loss.shape, (2,)) @@ -54,12 +60,13 @@ def test_random_loc_percentage(self): """Test clutter loss calculation with random location percentage.""" frequency = np.array([4000]) # MHz distance = np.array([1000]) # meters - + clutter_type = 'one_end' loss = self.clutter_loss.get_loss( frequency=frequency, distance=distance, loc_percentage="RANDOM", station_type=StationType.IMT_UE, + clutter_type=clutter_type ) self.assertTrue(0 <= loss <= 100) diff --git a/tests/test_sharc_geom.py b/tests/test_sharc_geom.py index 1c7e409ec..d0aee89dd 100644 --- a/tests/test_sharc_geom.py +++ b/tests/test_sharc_geom.py @@ -2,10 +2,12 @@ import numpy as np import numpy.testing as npt import shapely as shp +import shapely.vectorized from pathlib import Path from sharc.support.sharc_geom import generate_grid_in_multipolygon from sharc.support.sharc_utils import load_gdf +from sharc.satellite.utils.sat_utils import haversine class TestSharcGeom(unittest.TestCase): @@ -21,6 +23,7 @@ def test_generate_grid(self): """Test the generate_grid_in_multipolygon function.""" # approx a square # good only for small (lon, lat) values + rng = np.random.RandomState(seed=0xCAFEBABE) mx = 0.001 border = np.linspace(-mx / 2, mx / 2, 100) poly = shp.Polygon( @@ -30,16 +33,9 @@ def test_generate_grid(self): [(x, -mx / 2) for x in border[::-1]] ) # when passing too large values, the grid will have no points - # hexagon radius = x means square has at least side - # l >= min(x * sin(60deg), x * sin(30deg) + x/2) - # l >= min(sqrt(3) * x / 2, x) - # l >= sqrt(3) * x / 2 grid = generate_grid_in_multipolygon( poly, 1.01 * mx * 111e3 * 2 / np.sqrt(3)) self.assertEqual(len(grid[0]), 0) - grid = generate_grid_in_multipolygon( - poly, 0.9 * mx * 111e3 * 2 / np.sqrt(3)) - self.assertEqual(len(grid[0]), 1) # doing some packing inside the square # based on @@ -48,8 +44,12 @@ def test_generate_grid(self): pol_A = mx * mx * 111e3 * 111e3 for hx_r in np.arange(100, 12801, 5e2): hx_A = 3 * np.sqrt(3) * hx_r**2 / 2 - grid = generate_grid_in_multipolygon(poly, hx_r) + grid = generate_grid_in_multipolygon(poly, hx_r, True, rng) npt.assert_allclose(len(grid[0]), pol_A / hx_A, rtol=0.05, atol=10) + lons, lats = grid + self.assertTrue( + shp.vectorized.contains(poly, lons, lats).all() + ) gdf = load_gdf( self.countries_shapefile, @@ -59,16 +59,20 @@ def test_generate_grid(self): # 10km hexagon radius hx_r = 10e3 - grid = generate_grid_in_multipolygon(poly, hx_r) + grid = generate_grid_in_multipolygon(poly, hx_r, True, rng) pol_A = 851e4 * 1e6 hx_A = 3 * np.sqrt(3) * hx_r**2 / 2 br_grid_len = len(grid[0]) - npt.assert_allclose(br_grid_len, pol_A / hx_A, rtol=1e-5, atol=190) + npt.assert_allclose(br_grid_len, pol_A / hx_A, rtol=0.007, atol=190) + lons, lats = grid + self.assertTrue( + shp.vectorized.contains(poly, lons, lats).all() + ) poly = gdf[gdf["NAME"] == "Chile"]["geometry"].values[0] - grid = generate_grid_in_multipolygon(poly, hx_r) + grid = generate_grid_in_multipolygon(poly, hx_r, True, rng) pol_A = 756626 * 1e6 hx_A = 3 * np.sqrt(3) * hx_r ** 2 / 2 @@ -76,12 +80,34 @@ def test_generate_grid(self): # causes noticeable distortion and differing results... cl_grid_len = len(grid[0]) npt.assert_allclose(cl_grid_len, pol_A / hx_A, rtol=0.08, atol=230) + lons, lats = grid + self.assertTrue( + shp.vectorized.contains(poly, lons, lats).all() + ) # generating grid for both Chile and Brazil at the same time poly = shp.ops.unary_union(gdf["geometry"].values) - grid = generate_grid_in_multipolygon(poly, hx_r) + grid = generate_grid_in_multipolygon(poly, hx_r, True, rng) + + npt.assert_allclose(len(grid[0]), cl_grid_len + br_grid_len, rtol=0.0011) + lons, lats = grid + self.assertTrue( + shp.vectorized.contains(poly, lons, lats).all() + ) + + from scipy.spatial import cKDTree + + lons, lats = grid + pts = np.column_stack([lons, lats]) # shape (N, 2) + tree = cKDTree(pts) + dists, idxs = tree.query(pts, k=2) + # idxs[:,0] are the points themselves + # idxs[:,1] are the nearest neighbor + near_idxs = idxs[:, 1] + nn_dists = haversine(lons, lats, lons[near_idxs], lats[near_idxs]) - self.assertEqual(len(grid[0]), cl_grid_len + br_grid_len) + expected_dist = np.sqrt((np.sqrt(3) / 2 * hx_r)**2 + (3 / 2 * hx_r)**2) + npt.assert_allclose(nn_dists, expected_dist, rtol=0.02) if __name__ == '__main__': diff --git a/tests/test_sharc_utils.py b/tests/test_sharc_utils.py new file mode 100644 index 000000000..f6a1ae15b --- /dev/null +++ b/tests/test_sharc_utils.py @@ -0,0 +1,65 @@ +import unittest + +from sharc.support.sharc_utils import clip_angle + + +class StationTest(unittest.TestCase): + """Unit tests for some utilities.""" + + def setUp(self): + """ + setup + """ + pass + + def test_clip_angle(self): + """ + Testing in range for non wrapping around angles + """ + for a in [90, 100, 180, -180, -135.01]: + self.assertEqual(clip_angle(a, 0, 90), 90) + for a in [-134.99, -90, -45, 0]: + self.assertEqual(clip_angle(a, 0, 90), 0) + for a in [0, 45, 90]: + self.assertEqual(clip_angle(a, 0, 90), a) + + for a in [-90, -80, 0, 44.99]: + self.assertEqual(clip_angle(a, -180, -90), -90) + for a in [45.01, 90, 135, 180, -180]: + self.assertEqual(clip_angle(a, -180, -90), -180) + for a in [-180, -135, -90]: + self.assertEqual(clip_angle(a, -180, -90), a) + + for a in [180, -180, -135, -90.01]: + self.assertEqual(clip_angle(a, 0, 180), 180) + for a in [-89.99, -45, 0]: + self.assertEqual(clip_angle(a, 0, 180), 0) + for a in [0, 45, 90, 135, 180]: + self.assertEqual(clip_angle(a, 0, 180), a) + + """ + Testing in range for wrapping around angles + """ + for a in [-90, -80, 0, 44.99]: + self.assertEqual(clip_angle(a, 180, -90), -90) + for a in [45.01, 90, 135, 180, 180]: + self.assertEqual(clip_angle(a, 180, -90), 180) + for a in [-180, -135, -90]: + self.assertEqual(clip_angle(a, 180, -90), a) + + for a in [-180, -135, -90.01]: + self.assertEqual(clip_angle(a, 0, -180), -180) + for a in [-89.99, -45, 0]: + self.assertEqual(clip_angle(a, 0, -180), 0) + for a in [0, 45, 90, 135, -180]: + self.assertEqual(clip_angle(a, 0, -180), a) + + self.assertEqual(clip_angle(91, 180, 0), 180) + self.assertEqual(clip_angle(89, 180, 0), 0) + + self.assertEqual(clip_angle(91, -180, 0), -180) + self.assertEqual(clip_angle(89, -180, 0), 0) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_simulation_uplink.py b/tests/test_simulation_uplink.py index 9bf681ecf..def02bb60 100644 --- a/tests/test_simulation_uplink.py +++ b/tests/test_simulation_uplink.py @@ -868,7 +868,7 @@ def test_beamforming_gains(self): # Physical pointing angles self.assertEqual(self.simulation.bs.antenna[0].azimuth, 0) self.assertEqual(self.simulation.bs.antenna[0].elevation, -10) - self.assertEqual(self.simulation.bs.antenna[1].azimuth, 180) + self.assertEqual(self.simulation.bs.antenna[1].azimuth, -180) self.assertEqual(self.simulation.bs.antenna[0].elevation, -10) # Change UE pointing diff --git a/tests/test_spectral_mask_mss.py b/tests/test_spectral_mask_mss.py index feef1a9ec..f19bff26a 100644 --- a/tests/test_spectral_mask_mss.py +++ b/tests/test_spectral_mask_mss.py @@ -21,10 +21,10 @@ def test_power_calc(self): p_tx_density = -30 # dBW / Hz freq = 2190 band = 1 - p_tx = p_tx_density + 10 * np.log10(band * 1e6) + p_tx = p_tx_density + 10 * np.log10(band * 1e6) + 30 # reference bw is 4khz below 15GHz center freq - p_tx_over_4khz = p_tx_density + 10 * np.log10(4e3) + p_tx_over_4khz = p_tx_density + 10 * np.log10(4e3) + 30 # dBm/MHz spurious_emissions = -30 @@ -45,7 +45,7 @@ def test_power_calc(self): F = (f_offset - band / 2) / band * 100 - should_eq[2 * i] = p_tx_over_4khz - 40 * np.log10(F / 50 + 1) + 30 + should_eq[2 * i] = p_tx_over_4khz - 40 * np.log10(F / 50 + 1) eq[2 * i] = mask.power_calc(freq + f_offset + 0.5 * 4e-3, 4e-3) should_eq[2 * i + 1] = should_eq[2 * i] diff --git a/tests/test_topology_imt_mss_dc.py b/tests/test_topology_imt_mss_dc.py index 2207684a0..142af61ac 100644 --- a/tests/test_topology_imt_mss_dc.py +++ b/tests/test_topology_imt_mss_dc.py @@ -68,8 +68,6 @@ def test_initialization(self): self.assertEqual( self.imt_mss_dc_topology.num_sectors, self.params.num_beams) - self.assertEqual(len(self.imt_mss_dc_topology.orbits), - len(self.params.orbits)) def test_calculate_coordinates(self): """Test calculation of coordinates for the topology.""" @@ -109,9 +107,9 @@ def test_calculate_coordinates(self): # by default, satellites should always point to nadir (earth center) ref_earth_center = StationManager(1) - ref_earth_center.x = self.earth_center_x - ref_earth_center.y = self.earth_center_y - ref_earth_center.z = self.earth_center_z + ref_earth_center.x = self.earth_center_x.flatten() + ref_earth_center.y = self.earth_center_y.flatten() + ref_earth_center.z = self.earth_center_z.flatten() ref_space_stations = StationManager( self.imt_mss_dc_topology.num_base_stations) diff --git a/tests/test_topology_single_base_station.py b/tests/test_topology_single_base_station.py index 9e2687f2e..a5c6d6fc8 100644 --- a/tests/test_topology_single_base_station.py +++ b/tests/test_topology_single_base_station.py @@ -85,7 +85,7 @@ def test_azimuth(self): """Test calculation of azimuth for the topology.""" cell_radius = 1500 num_clusters = 1 - topology = TopologySingleBaseStation(cell_radius, num_clusters) + topology = TopologySingleBaseStation(cell_radius, num_clusters, [0.0]) topology.calculate_coordinates() self.assertEqual(topology.azimuth, 0) @@ -95,6 +95,47 @@ def test_azimuth(self): topology.calculate_coordinates() npt.assert_equal(topology.azimuth, np.array([0, 180])) + cell_radius = 1000 + num_clusters = 1 + topology = TopologySingleBaseStation(cell_radius, num_clusters, "random") + topology.calculate_coordinates() + self.assertEqual(len(topology.azimuth), num_clusters) + self.assertTrue(np.all((topology.azimuth >= 0) & (topology.azimuth < 360))) + + # Test that azimuth values are within the valid range [0, 360) + cell_radius = 1000 + num_clusters = 2 + topology = TopologySingleBaseStation(cell_radius, num_clusters, "random") + topology.calculate_coordinates() + self.assertEqual(len(topology.azimuth), num_clusters) + self.assertTrue(np.all((topology.azimuth >= 0) & (topology.azimuth < 360))) + + # Test that providing an azimuth list of incorrect length raises ValueError. + cell_radius = 1000 + num_clusters = 2 + topology = TopologySingleBaseStation(cell_radius, num_clusters, [0.0]) + with self.assertRaises(ValueError): + topology.calculate_coordinates() + + # Test that providing an invalid azimuth type raises ValueError + with self.assertRaises(ValueError): + TopologySingleBaseStation(cell_radius, num_clusters, "what?") + + # Test azimuth type and values for different cluster configurations. + cell_radius = 1000 + num_clusters = 1 + topology = TopologySingleBaseStation(cell_radius, num_clusters, [45.0]) + topology.calculate_coordinates() + self.assertTrue(isinstance(topology.azimuth, np.ndarray)) + self.assertEqual(topology.azimuth[0], 45.0) + + cell_radius = 1000 + num_clusters = 2 + azimuths = [30.0, 150.0] + topology = TopologySingleBaseStation(cell_radius, num_clusters, azimuths) + topology.calculate_coordinates() + npt.assert_equal(topology.azimuth, np.array(azimuths)) + if __name__ == '__main__': unittest.main()