Spectrum management and compatibility studies with Python

We developed the pycraf Python package, which provides functions and procedures for various tasks related to spectrum-management compatibility studies. This includes an implementation of ITU-R Rec. P.452, which allows to calculate the path attenuation arising from the distance and terrain properties between an interferer and the victim service. A typical example would be the calculation of interference levels at a radio telescope produced from a radio broadcasting tower. Furthermore, pycraf provides functionality to calculate atmospheric attenuation as proposed in ITU-R Rec. P.676. Using the rich ecosystem of scientific Python libraries and our pycraf package, we performed a large number of compatibility studies. Here, we will highlight a recent case study, where we analysed the potential harm that the next-generation cell-phone standard 5G could bring to observations at a radio observatory. For this we implemented a Monte-Carlo simulation to deal with the quasi-statistical spatial distribution of base stations and user devices around the radio astronomy station.


Introduction
The electromagnetic spectrum, especially at radio frequencies, is a limited resource, which has tremendous economic and scientific value for a large variety of services. Applications involve mobile communication networks, radio and TV broadcasting, RADAR, and safety distress signals. In contrast, natural sciences such as radio astronomy, Earth sensing (which is important for climatology), or weather forecasting rely heavily on interference-free spectral bands to conduct measurements with their high-sensitivity receivers. Especially, the radio-astronomy service (RAS) developed re-ceivers, often cryogenically cooled, that allow us to detect even the faintest signals from outer space. Today it is possible to observe Micro-Jansky sources, with a Jansky (Jy) being a unit of the spectral power flux density, which was introduced to honour Karl. G. Jansky the "founder" of the radio-astronomy field. It is 1 Jy = 10 −26 W m −2 Hz −1 .
With the increased utilization of the radio spectrum, it became essential to regulate its usage. Allocations of frequency bands are nowadays made by the International Telecommunication Union Radiocommunication Sector (ITU-R) in unanimous decisions by delegates of the administrations of its member states. Every two to four years, ITU-R organizes the World Radiocommunication Conference (WRC) where the so-called Radio Regulations are updated, which contain rules and procedures, and spectrum allocation tables. It is clear that such decisions and allocations have to be made in accordance to the technical needs and constraints of the concerned services. Therefore, in many regional and international organizations and institutions compatibility studies are performed, which analyse if, and possibly under which restrictions, a service can be granted a new allocation.
In the following, we present a new Python library, called pycraf, that can be used to carry out many of the recurring tasks, which occur in compatibility studies. This includes the estimation of path propagation losses and atmospheric dampening effects, the determination of satellite positions (given orbital parameters), or calculating antenna gain patterns. After a short introduction into compatibility studies in Section 2, the more important pycraf features are briefly introduced in Section 3.1 along with basic usage examples. A real-world application is discussed in Section 4, where the potential impact of the next-generation 5G mobile communication networks, which may acquire a new allocation at 24 GHz at the upcoming WRC-19, on the RAS is studied. Conclusions are presented in Section 5. arXiv:1805.11434v1 [astro-ph.IM] 29 May 2018 2 Spectrum management and compatibility studies Many compatibility studies follow a very simple recipe. The first step is usually to calculate the radiated power (or the power flux density), P tx , at the transmitter (Tx). Then one has to determine the appropriate interference threshold for the detected power, P rx , at the victim receiver (Rx). The ratio P tx /P rx constitutes the minimum shielding between interfering Tx and victim Rx, also called "minimum coupling loss" (MCL). In the final step the product of antenna gains and the path propagation loss (or path attenuation), i.e., the loss of power on the path to the receiver (Rx) is determined and compared to the MCL. That way one may find the minimum distance between a single interferer and victim receiver, which is often called separation distance. Since the path propagation loss also depends on the terrain height profile between Tx and Rx, the necessary separation distance will vary for different bearings, leading to non-circular exclusion zones around a victim receiver. In the case of a widespread deployment of many possible interferers one may need to calculate the expected total ("aggregated") interference levels from all interferers within a certain range around the victim receiver. In most cases this will require a statistical simulation of possible deployments in a given terrain. However, simplified numerical convolutions of path loss laws with the deployment distribution functions are often cheaper to compute and can also yield good estimates of the size of zones where interferers ought to be excluded.

Available software
The various tasks that have to be solved in typical compatibility studies are often similar. Nevertheless, major effort is necessary to develop a single tool, which could be used in all sorts of different scenarios.
One such tool is SEAMCAT (ECC Report 252), which was developed by the European Communications Committee (ECC) to be used by spectrum managers in CEPT 1 countries. It lets the user choose from a large amount of system parameters, deployment scenarios, and path propagation models -all with appropriate distribution functions -and applies Monte-Carlo sampling to calculate the desired output distributions, e.g., of the power levels detected at a receiver. Unfortunately, as of today, SEAMCAT can only perform socalled generic studies, i.e., the terrain profile between a transmitter and receiver is either neglected or approximated in a very coarse manner. For RAS stations in mountainous terrain this may lead to unrealistic results as the diffraction at hill or 1 The European Conference of Postal and Telecommunications Administrations (CEPT) is an organization where policy makers and regulators from 48 countries across Europe collaborate to harmonize telecommunication, radio spectrum and postal regulations. The European Communications Office (ECO) is the Secretariat of the CEPT. mountain tops can significantly increase the path attenuation between victim and interferer(s).
PathProfile is a publicly available 2 software tool for the prediction of local field strengths developed by Mike Willis. It features a GUI and uses terrain data for the calculation of path loss. As such, it may serve as the core of compatibility calculations, but will require additional specific interfacing for multi-interferer scenarios or aggregation studies.
A variety of MATLAB-based solutions to various spectrum management problems exists and is used by regulatory administrations, industry, and also by radio astronomers in spectrum management. Some of them have found approval by the ECC (e.g. ECC Report 247), but many are not in the public domain and less well known and are used mainly by individuals.

The pycraf package
While there are some existing software tools to aid with compatibility studies, they are either not suited for all relevant tasks (PathProfile), not easy to use from within a programming language (SEAMCAT), or are based on proprietary software (MATLAB-based scripts). Therefore, we decided to provide a number of functions and procedures, which were developed in the framework of several radio compatibility studies involving the Effelsberg observatory. The observatory is located in the Eifel mountains in Western Germany and operates a 100-m radio telescope. The code comes in the form of a Python library (a so-called package) that can easily be installed and used even by relatively inexperienced Python users. The package is named pycraf, open-source licensed (GPL-v3), and available on the Python package distribution server PyPI (Python Package Index 3 . The source code is furthermore hosted on the GitHub platform 4 , along with detailed documentation 5 , a bug tracker, and installation instructions.
We aim to make pycraf a community-driven project. Contributions to the source code or documentation are very welcome, but we are also happy to receive bug reports, feature requests and other proposals that would further improve the package.
The pycraf package makes use of the AstroPy package template providing a sophisticated build environment, which is based on setuptools, pytest, and sphinx. It allows the use of automated testing, makes incorporation of continuous integration services easy, and enables automated documentation generation. Because pycraf is available on the PyPI distribution server, one can use the Python tool pip to install the package into the local Python environment.
After successful installation, one can load and test the library with the following two statements: import pycraf pycraf.test (remote_data='any') With this remote_data option, pycraf will attempt to download terrain-height data from the NASA Shuttle Radar Topography Mission (SRTM Farr et al., 2007) for the calculation of path propagation loss in real terrain 6 .
In the following we will present an overview of the features included in pycraf and provide some basic examples for its use. A discussion of all aspects of pycraf and introduction to the various third-party libraries used in pycraf is beyond the scope of this paper and we refer to the online manual where a lot more information and a full API reference for each sub-package can be found. We also note that pycraf makes extensive use of the Astropy (physical) units module, which is well documented on the Astropy website 7 .

Conversions module
Most of the fundamental operations, e.g., converting the radiated total power into an equivalent power flux density at a given distance and assuming a certain antenna gain of the Tx, or calculating the total received power at an Rx from the power flux density, are covered by the conversions subpackage. Furthermore, several Decibel units are defined, e.g.: Listing 1. Conversions module.
The conversions sub-package offers a lot more and we refer to the online documentation for a detailed overview of the functionality.

Path attenuation
The largest part of the code base is currently dedicated to the implementation of the path propagation loss algorithm described in ITU-R Rec. P.452. As discussed in Section 2, the path attenuation calculation is one of the three fundamental steps in any compatibility study. Furthermore, it is the most critical aspect, as the losses can easily exceed 100 dB and typically uncertainties are large. The method in P.452 includes various effects that influence the path loss: 1. Line-of-sight (free-space) loss including correction terms for multipath and focussing effects, 2. Diffraction (at terrain features), 3. Tropospheric scatter, and 4. Anomalous propagation (ducting, reflection from elevated atmospheric layers).  (Farr et al., 2007).
Furthermore, a simple correction is provided to include additional losses caused by local obstacles ("clutter") at the end points of the propagation path. The details of the individual mechanisms are beyond the scope of this paper and we refer to the P.452 and references therein for further information.
We note, however, that the diffraction calculation was revised in version 15 of the P.452. Previously, it used the 3edge Deygout method, while the new algorithm is based on the so-called delta-Bullington method. Because this change can significantly affect the results of the diffraction loss in mountainous terrain (see Wilson and Salamon, 2008, and references therein), pycraf offers both methods, from the older version 14 and the current version 16, such that one can compare older study results with the predictions of the new method.
While most of pycraf is pure Python code, for performance reasons the pathprof sub-package contains a large amount of Cython 8 code (Behnel et al., 2011). Even then, the P.452 algorithm is computationally demanding and if one wants to create a large map of path attenuations around a transmitter, the calculation would be relatively slow. Allowing for an insignificant loss of numerical accuracy, we developed a fast alternative for this use case, which is presented in Section 3.2.1.
Before the path propagation loss can be calculated, a terrain height profile has to be constructed from a suitable data source. The pycraf package has built-in access to SRTM data servers, which allows to obtain so-called tiles (1 • × 1 • ) on demand: Listing 2. Path profile module.
1 import astropy.units as u 2 import pycraf.conversions as cnv 3 import pycraf.pathprof as pp 8 Cython is an extension to the Python language that allows the compilation of code into so-called C extensions, which can then be imported into Python and often provide a significant speed gain compared to Python-only implementations. The details of this are beyond the scope of this paper. More information can be found on the Cython webpage: http://cython.org/. Of course, it is also possible to provide SRTM tiles in a local directory to avoid excessive network traffic. The returned variables contain the geographical longitudes and latitudes of the path, the distances and back-bearings of each point on the path with respect to the transmitter (which are inferred via Vincenty's formulae for the Earth ellipsoid; Vincenty, 1975), as well as the heights (above the WGS84 mean sea level).
Here, the code to plot the data is omitted, but Fig. 1 displays the result. It is noted, however, that the relevant pycraf functions can query SRTM data automatically; users do not necessarily need to provide the height profile themselves. Calculating the propagation loss between two points is a two-step process. First, an auxiliary object, PathProp, is constructed that contains the geometrical, environmental, and terrain parameters necessary for the further calculations: Here, omega, is the fraction of the path that is over sea, time_percent is the percentage of time for which the propagation losses will be larger than the returned values (see P.452 for details), and h tg,rg are the antenna heights of the Tx and Rx over ground, respectively. Furthermore, the clutter zone type can be specified for both end points. The returned Python tuple contains 1. L bfsg , the free-space loss, 2. L bd , the basic transmission loss caused by diffraction, 3. L bs , the tropospheric scatter loss, 4. L ba , the ducting/layer reflection loss, 5. L b , the complete path propagation loss, 6. L b,corr , as L b but with clutter correction, and 7. L, as L b but with clutter and gain correction.
An important property of the total propagation loss (L b and L) is that it is not simply the product (or sum, if logarithmic representation is used) of the components. Again, we refer to P.452 for more details. To give the reader an impression on how the various quantities behave as a function of distance from the transmitter, we have calculated the components for each point in the terrain height profile ( Fig. 1) and display them in Fig. 2. In contrast to the two-point calculation above, the antenna gains were neglected, as the effective gain towards each point along the path differs somewhat due to the different elevation angles of the path at the two stations. That however requires accounting for the respective antenna patterns, which is supported by pycraf, but beyond the scope of this simple example.

A fast map-making algorithm
With the above approach it is straightforward to produce a path loss map, e.g., around a given transmitter, by simply calculating the propagation loss between each pixel in the map and the Tx location. As the P.452 algorithm is relatively complex, this would however require a substantial amount of computing time. In addition, the construction of the geodesics (the shortest paths between two points on the Earth ellipsoid) and SRTM height profile extraction/interpolation are also computationally demanding. Thus, two additional techniques were introduced to provide a significant speed-up. First, the map computation was moved into  (Farr et al., 2007).
the Cython part of the code (avoiding the use of slower Python for-loops), where it is furthermore possible to make use of OpenMP 9 for automatic parallelisation of the code to improve execution times on machines with many CPU cores. The second improvement is algorithmic in nature and aims to speed-up the geometric part of the computation: the geodesics calculation and height profile generation, which is preparatory work before the actual P.452-algorithm is employed. A path from the map centre to a pixel on the edge of the map crosses many other pixels in-between. We calculate all geodesics between map edge pixels and the centre and use a caching technique to identify the edge-geodesic which comes closest to the pixel centre for each interior pixel. This essentially transforms the problem from O(n 2 ) to O(n). Unfortunately, the P.452 algorithm itself can not be treated with a similar technique, because the path elevation angles (and thus also if a path is line-of-sight or trans-horizon) are unique for each pixel as they depend on the terrain height profile between the pixel and the map centre.
Using the attenuation map-making functionality is again a two-step process. First, a helper object is created that caches the terrain and path parameters, which is then fed into the P.452 algorithm: Listing 3. Fast path propagation loss map-making.
1 import astropy.units as u 2 import pycraf.conversions as cnv 3 import pycraf.pathprof as pp   (Farr et al., 2007). Again, we omit the lengthy plot code, but the result of this is visualized in Fig. 3. The hprof_cache object is a Python dict-like object containing the cached parameters and could also be edited by the user before calling atten_map_fast. A typical example of such editing would be the assignment of a different clutter zone to each of the pixels -given that the user has access to a suitable database containing the true clutter types.
The pathprof sub-package also provides a large amount of other utility functions, e.g., to plot terrain maps or export maps to a Geographic Information System (GIS) file format such as the Keyhole Markup Language (KML); see online documentation for further information.

Protection thresholds
Since pycraf was developed for compatibility studies involving radio astronomy, the current protection subpackage only offers thresholds relevant in that context. Most importantly, this covers the tables with protection thresholds given in ITU-R Rec. RA.769: Values from this table can then subsequently be used to compare received power levels of a potential interferer with the permitted power levels. Almost all compatibility studies that involve RAS refer to the RA.769 methodology or thresholds in one form or another.

Antenna patterns
In many realistic compatibility studies, the antenna patterns of transmitter and receiver have to be accounted for. This can be more complex than what might be naively assumed. For example, the base stations and partly also the user equipment of the next-generation mobile communications standard 5G will utilize phased array antennas. To calculate the effective gain for the propagation path, one must account for the direction of the formed beam, the orientation (and rotation) of the antenna frame normal vector, and the relative position of the propagation path with respect to the antenna frame. We will return to this in Section 4. In pycraf several (simplified) antenna models are provided that implement various ITU-R recommendations, which shall be used for spectrum management studies. As an example, the fairly complex 5G composite antenna patterns can be obtained in the following manner: Listing 5. 5G composite antenna patterns.
1 import numpy as np 2 import astropy.units as u 3 import pycraf.conversions as cnv 4 import pycraf.antenna as ant 5 6 azims = np.arange (-180, 180 where the azim_i and elev_j angles define the pointing of the formed beam in the antenna frame coordinates (azimuth and elevation). In Fig. 4

Atmospheric attenuation
Once the propagation path exceeds a certain length, e.g., for communication links to satellites or for radio-astronomical observation of extra-terrestrial sources, the attenuation by the gases in the atmosphere of the Earth becomes important, especially at higher frequencies. Again, there exists a ITU-R recommendation, P.676, which was implemented in the atm sub-package. It can be used to predict the attenuation as a function of frequency, path geometry, and atmospheric parameters. The two molecules that play a major role are oxygen (so-called dry component) and water (wet component), both introducing spectral-line emission and absorption. Furthermore, oxygen also produces the continuous nonresonant Debye-spectrum (below 10 GHz). At frequencies above 100 GHz continuous attenuation by nitrogen starts to play a role as well.
The calculation of atmospheric attenuation is performed via a ray-tracing approach. Let us assume a radio signal enters Earth's atmosphere and travels to a receiving station close to the ground. First, the atmosphere is divided into layers (several hundreds or more) and for each of these the attenuation of the incident signal is calculated, such that the signal, which leaves a layer is somewhat weaker than the signal that entered a layer. At the same time, the layer itself emits radio waves having the aforementioned spectral properties. Each layer has different physical properties, such as temperature, pressure, water content, and refractivity. These "height profiles" have to be known to calculate the specific attenuation per layer. For ITU-R related studies, it is widely accepted to work with mean profiles, which were averaged over a year. The "standard profile" is provided in ITU-R Rec. P.835 along with five more specialised profiles, associated with the geographic latitude and season ("high-latitude summer/winter", "mid-latitude summer/winter", and "low latitude"). For example: Listing 6. Atmospheric attenuation.
1 import numpy as np 2 import astropy.units as u The resulting profiles are displayed in Fig. 5. ITU-R Rec. P.676 contains tables with the resonance lines of oxygen and water, necessary to calculate the specific attenuation for one layer (or propagation paths close to the ground, which stay in a single layer). The function atten_specific_annex1 can be used for that: The resulting spectrum is shown in Fig. 6. Calculating such spectra for each layer allows us to perform the raytracing for a slant path through the full atmosphere (with the function atten_slant_annex1):  In Fig. 7 the total atmospheric attenuation of a slant path for Summer/Winter high-latitude standard profiles is presented, along the so-called equivalent-black-body temperature, T ebb , of the full atmosphere with major contributions by the atmospheric emission itself.

Satellite two-line elements
For studies involving a satellite or a satellite network it is often necessary to determine their position relative to an Earth station. Satellite orbits vary slowly in time as a consequence of the gravitational influence of planetary system  Figure 6. Specific attenuation for a layer with a total pressure of 1013 hPa, a temperature of 290 K, and a specific water density of 7.5 g cm −3 according to ITU-R Rec. P.676. bodies, most importantly the Moon, the Sun, and Jupiter. Furthermore, the orbits are subject to drag forces induced by the outer layers of the atmosphere, especially, in the case of low-earth orbit (LEO) satellites. Satellites are therefore constantly monitored by stations on Earth to update their orbital parameters. These are provided in a standardized format, the so-called two-line element (TLE) sets, which are usually updated once per day and can be retrieved from the Space Track 10 or CelesTrak 11 websites. As an example, the TLE of the International Space Station (ISS) for September, 20, 2008 looks like this: ISS (ZARYA) 1 25544U 98067A 0 8 2 6 4 . 5 1 7 8 2 5 2 8 −.00002182 00000−0 −11606−4 0 2927 2 25544 5 1 . 6 4 1 6 2 4 7 . 4 6 2 7 0006703 1 3 0 . 5 3 6 0 3 2 5 . 0 2 8 8 1 5 . 7 2 1 2 5 3 9 1 5 6 3 5 3 7 The meaning of the entries is described e.g. on Wikipedia 12 . The fields 7 and 8 (line 1, columns 19-20 and 21-32), i.e., 08 and 264.51782528 encode the epoch (year and day of year) at which the other orbital parameters are valid. From this point in time, one can then propagate the orbit.
To account for the aforementioned distortions of orbits, the calculations behind the propagation algorithm are rather complex and pycraf internally uses the Python package sgp4 for this. Provided with a TLE string, sgp4 can calculate the position of the satellite in the Geocentric reference frame for any point in time. However, if the time difference to the TLE epoch is significantly larger than a day, the result may be completely wrong.
For compatibility studies one usually needs the satellite position in the antenna reference frame of the ground station. Therefore, the satellite sub-package of pycraf, which is otherwise just a thin wrapper around sgp4, adds some routines for this coordinate conversion. One only needs to 10 https://www.space-track.org 11 http://celestrak.com/ 12 https://en.wikipedia.org/wiki/Two-line_element_set   As another example, Fig. 8 shows the position of the IRID-IUM network satellites, as seen from an observer at the Effelsberg observatory. IRIDIUM satellites are LEOs and move fast over the visible sky. To visualize this, the azimuth and elevation of the satellites were computed for a time range of five minutes.

Compatibility of a 5G cell phone network and RAS
With the pycraf package introduced in Section 3.1 it is very easy to do simple compatibility studies. Especially, for a point-to-point analysis, there is not much more to it than the determination of the emitted power and the path loss and finally the comparison of the received power with the compatibility threshold level. These steps could easily be put into a stand-alone software package, perhaps with a graphical user interface. A need for additional functionality may however arise sometimes, or one may be faced with more complex scenarios right from the beginning. Therefore, pycraf was designed as a programming library. Although this makes the simplest projects slightly more complicated (because one always need to write some Python code), it is well-suited to deal with more challenging studies. One of these more complicated cases is discussed in this section.

Context
The mobile telecommunication industry wants to utilize more bandwidth for their next-generation mobile communications standard, which is called 5G. At the last world radio communications conference (WRC), an agenda item (AI 1.13) was proposed for the next WRC in 2019. It is entitled "Work plan and proposed liaison to Task Group 5/1 on spectrum needs for the terrestrial component of IMT 13 in the frequency range between 24.25 GHz and 86 GHz". Although there is currently some tendency among vendors and administrations to favour the lower part of the frequency range mentioned in the agenda item (the 24.25-27.5 GHz sub-band) we are required to study all possible sub-bands without a predetermined conclusion for AI 1.13.
The Committee on Radio Astronomy Frequencies (CRAF) -an expert committee of European radio astronomers involved in spectrum management -is submitting input documents with compatibility studies to various regional and global working groups (ITU-R TG 5/1, ECC PT 1). These studies attempt to determine the compatibility of the nextgeneration 5G equipment with radio-astronomical observations. Depending on the designated 5G band, protected RAS bands would in some cases be located in the out-of-band (OOB) or spurious domain with respect to the IMT carrier frequency, in other cases we may face the in-band sharing scenario, where both, Tx and Rx, frequency ranges overlap. The latter case usually demands larger separation distances between interfering stations (the IMT equipment) and the victim station (the radio telescope) to ensure that radioastronomical observations do not suffer interference. Because of the additional attenuation of the emitted power a compatibility is more easily obtained in the OOB or spurious cases.
Here, we describe how the compatibility calculations for one of these studies were made, focussing on the 24.25-27.5 GHz frequency band. RAS has a primary allocation for 23.6-24 GHz, hence, this is an example for a OOB-domain scenario. For compatibility studies on the ITU level, one often analyses the "generic" case, where no terrain height (above the Earth ellipsoid) is assumed, so that the results do not depend on the terrain details of specific sites. Of course, RAS stations in mountainous regions may have some level of natural shield from surrounding hill tops, so that genericstudy results are only indicative. Still, many RAS stations are located in regions with relatively flat terrain, e.g., the Westerbork Synthesis Radio telescope in the Netherlands. We emphasize that the results presented here sensitively depend on the chosen transmitter properties and deployment densities. The assumptions and values, which are used in the following, are preliminary, as some of relevant (ITU-R) recommendations and guidelines are not yet finalized. Therefore, one should be careful with drawing conclusions from our study.

Setup I: Transmitter and receiver properties
5G networks will likely be "classic" mobile networks, consisting of base stations (BS) and user equipment (UE), although alternatives such as mesh-network based topologies seem also viable. Both, BS and UE feature phased-array antenna patterns, as described in Section 3.4. The transmitter properties (to be used for compatibility calculations) are summarized in ITU-R Doc. TG 5/1-36 (Attachment 2). For BS, 8 × 8 individual antenna elements, each having a power output of 10 dBm per carrier (before Ohmic losses), are foreseen. UE will use smaller arrays, of 4 × 4 elements, but with the same power levels. For the spurious domain, a spectrum block edge mask of −38 dBc per carrier bandwidth (200 MHz) for BS and −32 dBc/carrier for UE is in place, such that the spectral power density within the RAS band is −13 dBm/MHz in both cases. For UE an additional socalled body loss of 4 dB has to be accounted for, because in many situations a person holding the user device (e.g., a smart phone) will provide additional shielding. Integrating the spectral power over the RAS band yields a total of 13 dBm. ITU-R Doc. TG 5/1-36 (Attachment 2) makes a distinction for different zone types, urban, suburban and suburban open-space regions with different number densities and antenna heights. BS antenna heights will be 15 m for suburban open-space zones and 6 m for the rest, while UE is always at 1.5 m height. Even though the 5G phased array antennas allow dynamic beam forming, ITU-R Rec. M.2101 advises to use the single-element patterns for studies in the spurious domain. We follow this advice here, although there are diverging opinions about it as some consider it unlikely that the signal phasing should be completely inefficient close to the carrier frequency. Fig. 9 shows the resulting antenna gain patterns for BS and UE, respectively. The composite antenna patterns will still be required (see Section 4.6), because they influence the link budget (or coupling loss) between BS and UE, which in turn determines the level of power control in the user devices -UE can increase or decrease the transmitted power based on the coupling loss for more efficient battery use.
On the receiver side, an isotropic antenna (0 dBi) is usually assumed for these kinds of studies, with a height of 50 m above ground, which would be typical for a 100-m class radio telescope. ITU-R Rec. RA.769 permits an interference power level of −165 dBm (in continuum mode). Therefore, the minimal coupling loss (MCL), i.e., the attenuation induced by the path propagation loss plus the gain or attenuation due to the Tx antenna gain, must be at least 178 dB. The MCL ultimately determines the minimal distance (also called separation distance) that a device must have from the Rx to avoid problems with RAS observations.

Path propagation loss
With the functionality contained in pycraf, it is very easy to determine the path loss as a function of distance between Tx and Rx. Here, a distinction is to be made between the sub-/urban and suburban open-space areas, because the antenna height plays a role for the propagation. The default clutter model of ITU-R Rec. P.452 is not be applied in studies related to agenda item AI 1.13, instead the use of the model provided in ITU-R Rec. P.2108 is mandatory. It depends only on frequency, distance, and a quantity called location percentage, p L , which is the percentage of all emitters producing the lowest clutter loss. For example, if p L = 2%, the returned value, L clutter , indicates that for 2% of all devices the actual clutter loss will fall short of L clutter . In other words, P.2108 provides the cumulative distribution function for the clutter. We also note that the distance dependence in ITU-R Rec. P.2108 is only relevant for very small values below ∼2 km.
In Figs. 10 and 11 the resulting path propagation loss for BS and UE is displayed. The rapid increase of the path loss between about 35 and 50 km is caused by the geometry of the path. There, the line-of-sight path (small distances) turns into a trans-horizon path, meaning that diffraction kicks in.

Single-interferer scenario
The single-interferer scenario is an extreme case, where only one IMT device is considered having maximally possible antenna gain (i.e., the antenna shall be directed towards the RAS station 14 ) and very low p L (which means the clutter loss will be unusually small). In reality, such a situation will rarely occur, but it is still valuable to study this case, because it already gives a rough indication about the magnitude of the necessary separation distances. Very often, the single- 14 In case of BS, the array antennas are installed with a down-tilt of −10 • , so that the single-element pattern maximum is usually not pointing towards the RAS station. Total Total w. clutter Total w. clutter 2% Figure 11. As Fig. 10 for UE.
interferer separation distances define a zone, where operation of the interfering service will be critical. Figures 12 and 13 contain the resulting power levels received from a single BS or UE with the isotropic RAS antenna. The RA.769 threshold is marked by a dashed red line. The crossing point of the curves with the threshold defines the necessary separation distance. We note that for IMT networks of the previous generations, a spurious power level of −30 dBm/MHz was in place in Europe. To allow for comparison, the figures also contain the results for this lower level.

Setup II: Deployment of 5G devices
More important than the single-interferer analysis is the calculation of the received aggregated power levels, which stem from the full IMT network as the sum of powers from each individual device. Here one aims to use a realistic scenario for the deployment of the IMT equipment to the best of one's knowledge. This not only involves modelling how the BS and UE are distributed in a given area, but one must also account for the quasi-random orientations that the BS and especially the user devices may have. The latter can e.g. be freely rotated with respect to the BS.
Both, BS and UE will use phased-array technology to form beams into the direction of the associated communication partner. Relatively high gains can be achieved this way, even with the small antenna apertures. BS, which are associated to several user devices simultaneously, use the time-division duplex (TDD) method to serve the UE. In each time slot, the phased-array weights are updated to form a beam to the currently active device. The BS antenna frames themselves are fixed, in contrast to the user devices, where the system needs to form the beam depending on the location relative to the BS and its current rotation state. The beam forming also has an impact on the power that is transmitted into the direction of the RAS station, because the beam direction will significantly change the side-lobe pattern of the antennas (compare Fig. 4). Due to the TDD, each user device is only active during 20% of the time. For BS the TDD factor is set to 80%. Furthermore, in reality the network as a whole never runs at full capacity. For scenarios such as the one analysed here, ITU-R Doc. TG 5/1-36 (Attachment 2) recommends to work with a network loading factor of 20%.
In ITU-R Doc. TG 5/1-36 (Attachment 2) typical number densities are provided, which shall be used for generic compatibility studies involving large areas. The densities differ for BS and UE, but also for the various zone types. The parametrization is done in terms of the percentage of the considered area that has housing, R b = 5%, and the fraction within this area that is expected to have 5G installations, The final ingredient is the location of the UE relative to the associated BS. ITU-R Doc. TG 5/1-92 (Annex 1) defines distribution functions for the distance (Rayleigh for urban/suburban; Log-Normal for suburban open-space) and azimuthal position (Normal, N (0 • , 30 • )) relative to the antenna normal vector 15 . The two distance distributions are displayed in Fig. 14. It should be noted that each BS features only one sector antenna, in accordance to ITU-R Doc. TG 5/1-36 (Attachment 2).
In ITU-R Rec. M.2101 several possible deployment topologies are discussed, e.g., hexagonal or Manhattan-style grid layouts. In the particular case that is analysed here, the network topology can be neglected because one needs to average over a very large region, such that the aggregated power at the RAS station will be dominated by the (constant) 15 The azimuthal distribution is clipped for angles larger than deployment densities defined in ITU-R Doc. TG 5/1-36 (Attachment 2).

Aggregated power levels
We use a Monte-Carlo sampling approach to calculate the aggregated powers from the IMT network. This means that random values for BS and UE positions and rotations are generated from the distribution functions discussed in Section 4.5. This process is repeated many times and for each realization the total power received at the RAS station is calculated. From the ensemble of received powers (one value per trial run), one can then study distribution properties such as the median or average power level, but also the uncertainties (spread about the median).
The first step in our aggregation analysis is to generate a map containing zone types. A sufficiently large box of size 400 km × 400 km, with grid cells of size 0.5 km × 0.5 km is created. For each cell, a random number from a uniform distribution, U(0, 1), is drawn. Cells with values between 1 − (R urb a + R sub a )R b and 1 − R urb a R b are classified as suburban, while values larger than 1 − R urb a R b define urban areas. (Regions with housings are defined by values exceeding 1 − R b .) Figure 15 shows a zoom-in to the central part of the simulated box, with the random zone types assigned to it.
One could argue that such uniform-density deployment is not very realistic. Therefore, as an alternative we develop an algorithm that produces zones, which better resemble a real-world distribution of urban and suburban areas (inspired by the situation around the Effelsberg observatory). For this, some kind of clustering has to be introduced. One possibility to achieve that, is by smoothing the random numbers with one or more Gaussian filter kernels before assigning the zone types. The smoothing introduces a correlation be-tween neighbouring pixels. The larger the smoothing radius, the stronger the clustering of zone types will be. We find that a combination of three different kernels with scales σ k = 2, 5, and 15 km and relative amplitudes of 30%, 30%, and 40% can produce a zone-type map, which seems to be quite realistic; see Fig. 16.
Into the grid cells, BS are now sampled according to the zone types and given number densities. For each BS, a random azimuthal orientation of the antenna frame is assigned, and the tilting angle in elevation is considered. UE has to be handled differently, because several user devices can be associated to each BS, and the locations of the UE must be in the respective forward cone of the BS. Therefore, for each BS, user devices are sampled into the forward cone, until the total number of UE fits the desired global number density (per zone). UE can be freely rotated, with the only restriction, that the angular distance between antenna normal vector and direction to the associated BS shall not exceed 60 • ; see ITU-R Doc. TG 5/1-92 (Annex 1). Such a general rotation can for example be constructed from three Euler-angle rotations. Initially, the UE antenna normal shall point to the BS. We first rotate randomly about the UE-BS vector, then about the axis that is perpendicular to the UE-BS vector, is parallel to the x-y plane (the ground) and includes the UE position. Finally, another rotation about the UE-BS vector is applied. The first and last rotation angles are drawn uniformly from the interval [−180 • , 180 • ], while the second rotation has random angles from [−60 • , 60 • ]. The resulting rotation has the desired properties.
The next step is to calculate the effective antenna gains of each 5G device towards the RAS station. Since the positions and orientations of all devices are known, one can easily determine the azimuth and elevation of the UE-BS vectors, as well as of the direction to the RAS Rx, in the (rotated) antenna frames. Strictly speaking, if the singleelement antenna patterns are used, the UE-BS vectors are irrelevant. However, they will be needed for the power control algorithm, which is discussed below. Using the antenna gain formula provided in ITU-R Rec. M.2101 (implemented in pycraf.antenna), the effective gains follow directly. For BS, there is not a single resulting gain value, but several -one for each user device in the forward cone. Therefore, the average value is calculated and applied in the following. Figure 17 visualizes the geometry and the determined effective gains for a small sub-region that includes three BS and their associated UE. The red lines mark which user device belongs to which BS. The black arrows are the antenna normal vectors (projected into the x-y plane, i.e., shorter arrows have larger z components), while grey arrows point to the RAS station. The more the black and grey arrows are aligned, the higher the effective gain. If the composite patterns would be used, the arrows additionally had to be aligned with the red BS-UE vectors to produce highest effective gain. However, the beam forming makes the latter alignment less important. Only if the beam direction has large angular sepa- ration from the antenna normal, the achievable gain towards the beam direction drops significantly. However, the specific side-lobe patterns play an important role in the compositeantenna case. The path propagation loss to the RAS Rx is then computed for each 5G device, as discussed in Section 4.3, such that all necessary quantities to determine the total power level are now known. There is one last aspect, which makes the calculation slightly more complex: the power control algorithm mentioned in Section 4.2. For this, the link budget between BS and UE has to be determined (according to 3GPP TR 38.901, UMi -Street Canyon scenario). The details are beyond the scope of this paper. If the coupling loss between BS and UE is low, a user device may decrease the transmitter power to save battery power (according to a formula provided in ITU-R Rec. M.2101). The power control leads to a somewhat lower aggregated UE power, but since BS have a much larger impact on the total power (and are not subject to power control), the change is insignificant.
The aggregated powers received at the RAS station exceed the ITU-R Rec. RA.769 thresholds in every single Monte-Carlo run. This is not unusual because there is always a high chance that some devices get placed too close to the RAS Rx. To find out, how big an exclusion zone (i.e., the lower limit of separation distances to each device) would need to be, we show cumulative plots of the aggregated powers in Figs. 18 and 19. For increasing radii, r i , devices within a sphere with this radius are left out from the aggregation. The larger r i , the lower becomes the received power, as only more distant devices contribute. In total, one thousand runs were carried out. For each of these realizations, Fig. 18 contains the cumulative distribution for both, BS and UE, in the uniform-density scenario. Figure 19 shows the results for the clustered-density approach, revealing much higher scatter.
In order to determine a single value for the necessary exclusion zone radius, we compute distribution percentiles at 50% (the median) and 98% (the top 2% level); see Figs. 20 and 21. As can be seen, BS produce on average more than 10 dB higher aggregation levels than UE. Furthermore, there is a steep regime between 30 km and 60 km, where the power level drops quickly. This behaviour can be attributed to the path propagation loss, which shows a similar kink (compare Figs. 10 and 11). Hence, the transmitted power levels could be increased substantially, without the exclusion zone size having to grow quickly (unless one would not want to accept separation distances of ∼60 km in the first place).

Conclusions
The pycraf package proves to be a valuable and versatile tool to assist spectrum managers with accomplishing com-