Repository: bwinkel/pycraf (language: python)

Issue:
Matplotlib 3.3.3 -- 03c_attenuation_maps failing
When running example notebook 03c_attenuation_maps, the color bar fails to render correctly with the proper tick placement, labels, and number of ticks. This occurs in code cell 8:
`ValueError: The number of FixedLocator locations (5), usually from a call to set_ticks, does not match the number of ticklabels (7).`

I am using Matplotlib 3.3.3 on Ubuntu 18 in python 3.6.9. 


Discussion / hints:
The notebooks work under matplotlib 3.2.2
Thanks a lot for the hint. New matplotlib versions keep being maintenance heavy... I'll fix this in the next release.
I actually found the fix for this. It was painfully simple. The default vmax for the cim object was 950, but the colorbar was setting ticks up to 1350, so changing those to match (either one) fixed the problem.

That being said, it appears the pathprof wrapper function terrain_cmap_factory only takes an optional vmax argument but no vmin. This now causes a warning from matplotlib (attached). When using the color norm created for the terrain, matplotlib now wants vmin and vmax to be set when creating the norm.
![Screen Shot 2021-07-30 at 1 02 17 PM](https://user-images.githubusercontent.com/29465998/127687407-7523cee1-2b12-4cd3-b206-e9721a67232f.png)

Good to know. I stumbled upon the first problem myself recently, when I updated matplotlib. I'll take care of it asap.

Fix patch:
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 7e27f0ae8..f89c8cbf7 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -107,11 +107,11 @@ jobs:
       #   numpy.version: '1.13.1'
       #   astropy.version: '2'
       #   sphinx.version: '<3.0'
-      mac36:
-        imageName: 'macos-10.14'
-        python.version: '3.6'
-        numpy.version: '1.13.1'
-        astropy.version: '3'
+      # mac36:
+      #   imageName: 'macos-10.14'
+      #   python.version: '3.6'
+      #   numpy.version: '1.13.1'
+      #   astropy.version: '3'
       mac37:
         imageName: 'macos-10.14'
         python.version: '3.7'
diff --git a/docs/antenna/index.rst b/docs/antenna/index.rst
index 7f86e0375..99e017ea6 100644
--- a/docs/antenna/index.rst
+++ b/docs/antenna/index.rst
@@ -19,7 +19,7 @@ Using `pycraf.antenna`
 Antenna patterns for fixed and mobile service (ITU-R Rec. F.1336)
 -----------------------------------------------------------------
 For fixed link and mobile communications (IMT), several antenna patterns are
-provided in `ITU-R Rec F.1336-4 <https://www.itu.int/rec/R-REC-F.1336-4-201402-I/en>`_.
+provided in `ITU-R Rec F.1336-5 <https://www.itu.int/rec/R-REC-F.1336-5-201901-I/en>`_.
 
 Below, you can find a typical use case:
 
@@ -69,9 +69,10 @@ tilt angles.
 
 .. note::
 
-    At the moment only the sectoral peak-side lobe pattern for the frequency
-    range between 400 MHz and 6 GHz is provided. This is the one, which can
-    be used for IMT-advanced (LTE) basestations. Other patterns will follow.
+    At the moment only the sectoral peak- and average side lobe pattern for
+    the frequency range between 400 MHz and 6 GHz are provided. This is the
+    one, which can be used for IMT-advanced (LTE) basestations. Other
+    patterns will follow.
 
 WRC-19 Agenda Item 1.13, IMT2020
 --------------------------------
diff --git a/docs/geospatial/index.rst b/docs/geospatial/index.rst
index 4c7915304..ebdd50721 100644
--- a/docs/geospatial/index.rst
+++ b/docs/geospatial/index.rst
@@ -120,6 +120,7 @@ output). One can specifically ask for the original units by doing::
 
     >>> import pycraf.geospatial as geo
     >>> import astropy.units as u
+    >>> from astropy.units.imperial import ft
     >>> import pyproj
 
     >>> proj_wgs84 = pyproj.Proj('+init=epsg:4326')
@@ -146,11 +147,11 @@ because we can use the `~astropy.units` conversion::
     >>> x, y  # doctest: +FLOAT_CMP
     (<Quantity 925806.5486332772 m>, <Quantity 216168.1432314818 m>)
 
-    >>> x.to(u.imperial.ft), y.to(u.imperial.ft)  # doctest: +FLOAT_CMP
+    >>> x.to(ft), y.to(ft)  # doctest: +FLOAT_CMP
     (<Quantity 3037416.9849743457 ft>, <Quantity 709211.6499186204 ft>)
 
     >>> transform = geo.transform_factory(3452, 4326)
-    >>> transform(3037416.985 * u.imperial.ft, 709211.650 * u.imperial.ft)  # doctest: +FLOAT_CMP
+    >>> transform(3037416.985 * ft, 709211.650 * ft)  # doctest: +FLOAT_CMP
     (<Quantity -92.10581908573734 deg>, <Quantity 30.447921938477027 deg>)
 
 Unfortunately, there seems to be no way to ask `pyproj` about which
diff --git a/notebooks/03a_path_propagation_basic.ipynb b/notebooks/03a_path_propagation_basic.ipynb
index abb8a5748..91da84387 100644
--- a/notebooks/03a_path_propagation_basic.ipynb
+++ b/notebooks/03a_path_propagation_basic.ipynb
@@ -241,7 +241,7 @@
     "    _heightmap,\n",
     "    origin='lower', interpolation='nearest',\n",
     "    cmap=terrain_cmap, norm=terrain_norm,\n",
-    "    vmin=vmin, vmax=vmax,\n",
+    "    # vmin=vmin, vmax=vmax,  # deprecated in newer matplotlib versions\n",
     "    extent=(_lons[0], _lons[-1], _lats[0], _lats[-1]),\n",
     "    )\n",
     "cbar = fig.colorbar(\n",
@@ -257,6 +257,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 1150, 50)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} m'.format, ctics), color='k')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
diff --git a/notebooks/03c_attenuation_maps.ipynb b/notebooks/03c_attenuation_maps.ipynb
index c8b473998..879d185ff 100644
--- a/notebooks/03c_attenuation_maps.ipynb
+++ b/notebooks/03c_attenuation_maps.ipynb
@@ -232,7 +232,7 @@
     "    _heightmap,\n",
     "    origin='lower', interpolation='nearest',\n",
     "    cmap=terrain_cmap, norm=terrain_norm,\n",
-    "    vmin=vmin, vmax=vmax,\n",
+    "    # vmin=vmin, vmax=vmax,  # deprecated in newer matplotlib versions\n",
     "    extent=(_lons[0], _lons[-1], _lats[0], _lats[-1]),\n",
     "    )\n",
     "cbar = fig.colorbar(\n",
@@ -248,6 +248,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 1350, 200)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} m'.format, ctics), color='k')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
@@ -431,6 +432,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 300, 20)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
@@ -563,6 +565,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 300, 20)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
diff --git a/notebooks/03f_working_with_clutter.ipynb b/notebooks/03f_working_with_clutter.ipynb
index 1460b3f49..358e7f151 100755
--- a/notebooks/03f_working_with_clutter.ipynb
+++ b/notebooks/03f_working_with_clutter.ipynb
@@ -1450,6 +1450,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 300, 20)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
@@ -1539,6 +1540,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 300, 20)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
@@ -1599,6 +1601,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 26, 5)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
diff --git a/notebooks/03g_pre+postprocessing.ipynb b/notebooks/03g_pre+postprocessing.ipynb
index 3c3d4cb1c..741bd24df 100644
--- a/notebooks/03g_pre+postprocessing.ipynb
+++ b/notebooks/03g_pre+postprocessing.ipynb
@@ -192,6 +192,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(90, 131, 20)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} m'.format, ctics), color='k')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
@@ -735,6 +736,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(vmin + 10, vmax, 20)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
@@ -1167,6 +1169,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(vmin + 10, vmax, 20)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
@@ -1236,6 +1239,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(vmin + 10, vmax, 20)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
diff --git a/notebooks/A01_case_study_fixed_link.ipynb b/notebooks/A01_case_study_fixed_link.ipynb
index f248ec71e..f23382e1f 100644
--- a/notebooks/A01_case_study_fixed_link.ipynb
+++ b/notebooks/A01_case_study_fixed_link.ipynb
@@ -338,7 +338,7 @@
     "    _heightmap,\n",
     "    origin='lower', interpolation='nearest',\n",
     "    cmap=terrain_cmap, norm=terrain_norm,\n",
-    "    vmin=vmin, vmax=vmax,\n",
+    "    # vmin=vmin, vmax=vmax,  # deprecated in newer matplotlib versions\n",
     "    extent=(_lons[0], _lons[-1], _lats[0], _lats[-1]),\n",
     "    )\n",
     "cbar = fig.colorbar(\n",
@@ -354,6 +354,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 1150, 200)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} m'.format, ctics), color='k')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
diff --git a/notebooks/B02_wind_turbine_terrain.ipynb b/notebooks/B02_wind_turbine_terrain.ipynb
index 2f9e1c73c..2d89423bb 100644
--- a/notebooks/B02_wind_turbine_terrain.ipynb
+++ b/notebooks/B02_wind_turbine_terrain.ipynb
@@ -197,6 +197,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 300, 20)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
@@ -257,6 +258,7 @@
     "    t.label2.set_visible(True)\n",
     "\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "clabels = ['{:.0f} dB'.format(ct) for ct in ctics]\n",
     "clabels[0] = '                             <-20 dB'  # OMG, this is ugly\n",
     "clabels[-1] = '>20 dB                             '\n",
diff --git a/notebooks/B03_wind_turbine_reflection.ipynb b/notebooks/B03_wind_turbine_reflection.ipynb
index 021d44e12..d05b7dce7 100644
--- a/notebooks/B03_wind_turbine_reflection.ipynb
+++ b/notebooks/B03_wind_turbine_reflection.ipynb
@@ -475,7 +475,7 @@
     "    _heightmap,\n",
     "    origin='lower', interpolation='nearest',\n",
     "    cmap=terrain_cmap, norm=terrain_norm,\n",
-    "    vmin=vmin, vmax=vmax,\n",
+    "    # vmin=vmin, vmax=vmax,  # deprecated in newer matplotlib versions\n",
     "    extent=(_lons[0], _lons[-1], _lats[0], _lats[-1]),\n",
     "    )\n",
     "cbar = fig.colorbar(\n",
@@ -491,6 +491,7 @@
     "    t.label2.set_visible(True)\n",
     "ctics = np.arange(0, 1150, 200)\n",
     "cbar.set_ticks(ctics)\n",
+    "ctics = cbar.get_ticks()\n",
     "cbar.ax.set_xticklabels(map('{:.0f} m'.format, ctics), color='k')\n",
     "ax.set_xlabel('Longitude [deg]')\n",
     "ax.set_ylabel('Latitude [deg]')\n",
diff --git a/pycraf/antenna/cyantenna.pyx b/pycraf/antenna/cyantenna.pyx
index 18e56a985..37d5fda9b 100644
--- a/pycraf/antenna/cyantenna.pyx
+++ b/pycraf/antenna/cyantenna.pyx
@@ -381,7 +381,7 @@ cdef float64_t  _G_hr(
         return G
 
 
-cdef float64_t  _G_vr(
+cdef float64_t  _G_vr_peak(
         float64_t x_v, float64_t k_v, float64_t k_p,
         float64_t theta_3db, float64_t G180
         ) nogil:
@@ -408,6 +408,33 @@ cdef float64_t  _G_vr(
         return G180
 
 
+cdef float64_t  _G_vr_avg(
+        float64_t x_v, float64_t k_v, float64_t k_a,
+        float64_t theta_3db, float64_t G180
+        ) nogil:
+
+    cdef:
+
+        float64_t x_k = sqrt(1.33 - 0.33 * k_v)
+        float64_t C = (
+            10 * log10(
+                (180. / theta_3db) ** 1.5 * (4 ** -1.5 + k_v) / (1 + 8 * k_a)
+                ) /
+            log10(22.5 / theta_3db)
+            )
+        float64_t lambda_kv = 12 - C * log10(4.) - 10 * log10(4 ** -1.5 + k_v)
+
+    if x_v < x_k:
+        return -12 * x_v ** 2
+    elif x_k <= x_v and x_v < 4:
+        return -15 + 10 * log10(x_v ** -1.5 + k_v)
+    elif 4 <= x_v and x_v < 90 / theta_3db:
+        return -lambda_kv - 3 - C * log10(x_v)
+    else:
+        # x_v == 90 / theta_3db
+        return G180
+
+
 cdef float64_t _imt_advanced_sectoral_peak_sidelobe_pattern(
         float64_t azim, float64_t elev,
         float64_t G0, float64_t phi_3db, float64_t theta_3db,
@@ -477,7 +504,7 @@ cdef float64_t _imt_advanced_sectoral_peak_sidelobe_pattern(
     return (
         G0 +
         _G_hr(x_h, k_h, G180) +
-        R * _G_vr(x_v, k_v, k_p, theta_3db, G180)
+        R * _G_vr_peak(x_v, k_v, k_p, theta_3db, G180)
         )
 
 
@@ -557,6 +584,155 @@ def imt_advanced_sectoral_peak_sidelobe_pattern_cython(
     return it.operands[10]
 
 
+cdef float64_t _imt_advanced_sectoral_avg_sidelobe_pattern(
+        float64_t azim, float64_t elev,
+        float64_t G0, float64_t phi_3db, float64_t theta_3db,
+        float64_t k_a, float64_t k_h, float64_t k_v,
+        float64_t tilt_m, float64_t tilt_e,
+        ) nogil:
+
+    cdef:
+
+        float64_t beta = 0.
+        float64_t tmp = 0.
+
+        float64_t azim_rot = 0.
+        float64_t elev_rot = 0.
+
+        float64_t x_h = 0.
+        float64_t x_v = 0.
+
+        float64_t G180 = 0.
+        float64_t R = 0.
+
+    if tilt_m != 0.:
+
+        azim *= DEG2RAD
+        elev *= DEG2RAD
+        beta = tilt_m * DEG2RAD
+
+        elev_rot = asin(
+            sin(elev) * cos(beta) +
+            cos(elev) * cos(azim) * sin(beta)
+            )
+
+        tmp = (
+            -sin(elev) * sin(beta) +
+            cos(elev) * cos(azim) * cos(beta)
+            ) / cos(elev_rot)
+        if tmp > 1.:
+            tmp = 1.
+        if tmp < -1.:
+            tmp = -1.
+
+        azim_rot = acos(tmp)
+
+        azim = RAD2DEG * azim_rot
+        elev = RAD2DEG * elev_rot
+
+    if tilt_e != 0.:
+
+        tmp = elev + tilt_e
+        if tmp >= 0.:
+            elev_rot = 90 * tmp / (90 + tilt_e)
+        else:
+            elev_rot = 90 * tmp / (90 - tilt_e)
+
+        elev = elev_rot
+
+    G180 = -15 + 10 * log10(1 + 8 * k_a) - 15 * log10(180. / theta_3db)
+
+    x_h = fabs(azim) / phi_3db
+    x_v = fabs(elev) / theta_3db
+
+    R = (
+        (_G_hr(x_h, k_h, G180) - _G_hr(180. / phi_3db, k_h, G180)) /
+        (_G_hr(0., k_h, G180) - _G_hr(180. / phi_3db, k_h, G180))
+        )
+
+    return (
+        G0 +
+        _G_hr(x_h, k_h, G180) +
+        R * _G_vr_avg(x_v, k_v, k_a, theta_3db, G180)
+        )
+
+
+def imt_advanced_sectoral_avg_sidelobe_pattern_cython(
+        azim, elev,
+        G0, phi_3db, theta_3db,
+        k_a, k_h, k_v,
+        tilt_m, tilt_e,
+        gain=None,
+        ):
+    '''
+    Parallelized IMT advanced (LTE) antenna pattern (sectoral, avg side-lobe)
+    '''
+
+    cdef:
+
+        np.ndarray[float64_t] _azim, _elev
+        np.ndarray[float64_t] _G0, _phi_3db, _theta_3db
+        np.ndarray[float64_t] _k_a, _k_h, _k_v
+        np.ndarray[float64_t] _tilt_m, _tilt_e
+        np.ndarray[float64_t] _gain
+
+        int i, size
+
+    it = np.nditer(
+        [
+            azim, elev,
+            G0, phi_3db, theta_3db,
+            k_a, k_h, k_v,
+            tilt_m, tilt_e,
+            gain
+            ],
+        flags=['external_loop', 'buffered', 'delay_bufalloc'],
+        op_flags=[
+            ['readonly'], ['readonly'], ['readonly'], ['readonly'],
+            ['readonly'], ['readonly'], ['readonly'], ['readonly'],
+            ['readonly'], ['readonly'],
+            ['readwrite', 'allocate'],
+            ],
+        op_dtypes=[
+            FLOAT64, FLOAT64, FLOAT64, FLOAT64,
+            FLOAT64, FLOAT64, FLOAT64, FLOAT64,
+            FLOAT64, FLOAT64,
+            FLOAT64,
+            ]
+        )
+
+    # it would be better to use the context manager but
+    # "with it:" requires numpy >= 1.14
+
+    it.reset()
+
+    for itup in it:
+        _azim = itup[0]
+        _elev = itup[1]
+        _G0 = itup[2]
+        _phi_3db = itup[3]
+        _theta_3db = itup[4]
+        _k_a = itup[5]
+        _k_h = itup[6]
+        _k_v = itup[7]
+        _tilt_m = itup[8]
+        _tilt_e = itup[9]
+        _gain = itup[10]
+
+        size = _gain.shape[0]
+
+        for i in prange(size, nogil=True):
+
+            _gain[i] = _imt_advanced_sectoral_avg_sidelobe_pattern(
+                _azim[i], _elev[i],
+                _G0[i], _phi_3db[i], _theta_3db[i],
+                _k_a[i], _k_h[i], _k_v[i],
+                _tilt_m[i], _tilt_e[i],
+                )
+
+    return it.operands[10]
+
+
 # d_wlen = diameter / wavelength
 cdef float64_t _fl_pattern_2_1(
         float64_t phi, float64_t d_wlen, float64_t G_max
diff --git a/pycraf/antenna/imt.py b/pycraf/antenna/imt.py
index 45c972285..553cd384a 100644
--- a/pycraf/antenna/imt.py
+++ b/pycraf/antenna/imt.py
@@ -10,6 +10,7 @@
 from .cyantenna import imt2020_single_element_pattern_cython
 from .cyantenna import imt2020_composite_pattern_cython
 from .cyantenna import imt_advanced_sectoral_peak_sidelobe_pattern_cython
+from .cyantenna import imt_advanced_sectoral_avg_sidelobe_pattern_cython
 from .. import conversions as cnv
 from .. import utils
 
@@ -17,6 +18,7 @@
 __all__ = [
     'imt2020_single_element_pattern', 'imt2020_composite_pattern',
     'imt_advanced_sectoral_peak_sidelobe_pattern_400_to_6000_mhz',
+    'imt_advanced_sectoral_avg_sidelobe_pattern_400_to_6000_mhz',
     ]
 
 
@@ -334,7 +336,7 @@ def imt_advanced_sectoral_peak_sidelobe_pattern_400_to_6000_mhz(
     '''
     IMT advanced (LTE) antenna pattern (sectoral, peak side-lobe)
     for the frequency range 400 MHz to 6 GHz according to
-    `ITU-R Rec F.1336-4 <https://www.itu.int/rec/R-REC-F.1336-4-201402-I/en>`_
+    `ITU-R Rec F.1336-5 <https://www.itu.int/rec/R-REC-F.1336-5-201901-I/en>`_
     Section 3.1.1.
 
     Parameters
@@ -367,7 +369,7 @@ def imt_advanced_sectoral_peak_sidelobe_pattern_400_to_6000_mhz(
     Notes
     -----
     For typical values of :math:`k_p`, :math:`k_h`, and :math:`k_v` see
-    `ITU-R Rec F.1336-4 <https://www.itu.int/rec/R-REC-F.1336-4-201402-I/en>`_
+    `ITU-R Rec F.1336-5 <https://www.itu.int/rec/R-REC-F.1336-5-201901-I/en>`_
     Sections 3.1.1.1-3.
 
     For cases involving sectoral antennas with :math:`\\varphi_\\mathrm{3dB} \\lesssim 120^\\circ` the following
@@ -387,5 +389,80 @@ def imt_advanced_sectoral_peak_sidelobe_pattern_400_to_6000_mhz(
         )
 
 
+@utils.ranged_quantity_input(
+    azim=(-180, 180, apu.deg),
+    elev=(-90, 90, apu.deg),
+    G0=(None, None, cnv.dB),
+    phi_3db=(0, None, apu.deg),
+    theta_3db=(0, None, apu.deg),
+    k_a=(0, None, cnv.dimless),
+    k_h=(0, None, cnv.dimless),
+    k_v=(0, None, cnv.dimless),
+    tilt_m=(-90, 90, apu.deg),
+    tilt_e=(-90, 90, apu.deg),
+    strip_input_units=True, output_unit=cnv.dB
+    )
+def imt_advanced_sectoral_avg_sidelobe_pattern_400_to_6000_mhz(
+        azim, elev,
+        G0, phi_3db, theta_3db,
+        k_a, k_h, k_v,
+        tilt_m=0., tilt_e=0.,
+        ):
+    '''
+    IMT advanced (LTE) antenna pattern (sectoral, peak side-lobe)
+    for the frequency range 400 MHz to 6 GHz according to
+    `ITU-R Rec F.1336-5 <https://www.itu.int/rec/R-REC-F.1336-5-201901-I/en>`_
+    Section 3.1.2.
+
+    Parameters
+    ----------
+    azim, elev : `~astropy.units.Quantity`
+        Azimuth/Elevation [deg]
+    G0 : `~astropy.units.Quantity`
+        Antenna maximum gain [dBi]
+    phi_3db, theta_3db : float
+        3-dB beamwidth in the azimuth/elevation plane [degrees]
+    k_a : float
+        Parameter which accomplishes the relative minimum gain for
+        average side-lobe patterns [dimless]
+    k_h : float
+        Azimuth pattern adjustment factor based on leaked power
+        (:math:`0 \leq k_h \leq 1`) [dimless]
+    k_v : float
+        Elevation pattern adjustment factor based on leaked power
+        (:math:`0 \leq k_v \leq 1`) [dimless]
+    tilt_m : float
+        Mechanical tilt angle (downwards) [deg]
+    tilt_e : float
+        Electrical tilt angle (downwards) [deg]
+
+    Returns
+    -------
+    G : `~astropy.units.Quantity`
+        Antenna pattern [dB]
+
+    Notes
+    -----
+    For typical values of :math:`k_p`, :math:`k_h`, and :math:`k_v` see
+    `ITU-R Rec F.1336-4 <https://www.itu.int/rec/R-REC-F.1336-4-201402-I/en>`_
+    Sections 3.1.1.1-3.
+
+    For cases involving sectoral antennas with :math:`\\varphi_\\mathrm{3dB} \\lesssim 120^\\circ` the following
+    formula can be used to calculate :math:`\\vartheta_\\mathrm{3dB}`:
+
+    .. math::
+
+        \\vartheta_\\mathrm{3dB}=\\frac{31000\\times10^{-0.1G_0}}{\\varphi_\\mathrm{3dB}}
+
+    '''
+
+    return imt_advanced_sectoral_avg_sidelobe_pattern_cython(
+        azim, elev,
+        G0, phi_3db, theta_3db,
+        k_a, k_h, k_v,
+        tilt_m, tilt_e,
+        )
+
+
 if __name__ == '__main__':
     print('This not a standalone python program! Use as module.')
diff --git a/pycraf/pathprof/__init__.py b/pycraf/pathprof/__init__.py
index f32d532ba..b0bfaccac 100644
--- a/pycraf/pathprof/__init__.py
+++ b/pycraf/pathprof/__init__.py
@@ -14,6 +14,7 @@
 from .heightprofile import *
 from .helper import *
 from .propagation import *
+from .imt import *
 from .srtm import *
 
 _clutter_table = '''
diff --git a/pycraf/pathprof/cyimt.pyx b/pycraf/pathprof/cyimt.pyx
new file mode 100644
index 000000000..a7c2207b4
--- /dev/null
+++ b/pycraf/pathprof/cyimt.pyx
@@ -0,0 +1,362 @@
+#!python
+# -*- coding: utf-8 -*-
+# cython: language_level=3
+# cython: cdivision=True, boundscheck=False, wraparound=False
+# cython: embedsignature=True
+
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+cimport cython
+from cython.parallel import prange, parallel
+cimport numpy as np
+from numpy cimport PyArray_MultiIter_DATA as Py_Iter_DATA
+from libc.math cimport (
+    exp, sqrt, fabs, M_PI, sin, cos, tan, asin, acos, atan2, fmod, log, log10
+    )
+import numpy as np
+
+np.import_array()
+
+# __all__ = ['inverse', 'direct']
+
+
+cdef double NAN = np.NAN
+cdef double DEG2RAD = M_PI / 180.
+cdef double RAD2DEG = 180. / M_PI
+cdef double M_2PI = 2 * M_PI
+cdef double C = 3e8  # m/s
+
+
+cdef inline double _min(double a, double b) nogil:
+    if a < b:
+        return a
+    else:
+        return b
+
+
+cdef inline double _max(double a, double b) nogil:
+    if a > b:
+        return a
+    else:
+        return b
+
+
+cdef (double, double) _rural_macro_losses(
+        double fc_ghz, double d_2d,
+        double h_bs, double h_ue,
+        double W, double h,
+        ) nogil:
+    # all parameters (except freq) in meters;
+    # see 3GPP TR 38.901 version 16.1.0 Release 16;
+    # Table 7.4.1-1: Pathloss models; "RMa"
+
+    cdef:
+        double d_3d, d_bp, d_bp_3d
+        double PL1, PL2, PL_los, PL_nlos_prime, PL_nlos
+
+    d_3d = sqrt(d_2d ** 2 + (h_bs - h_ue) ** 2)
+    d_bp = M_2PI * h_bs * h_ue * fc_ghz * 1e9 / C
+    d_bp_3d = sqrt(d_bp ** 2 + (h_bs - h_ue) ** 2)
+
+    if d_2d < 10 or d_2d > 10000:
+        return (NAN, NAN)
+
+    if d_2d <= d_bp:
+        PL_los = PL1 = (
+            20 * log10(40 * M_PI * d_3d * fc_ghz / 3) +
+            _min(0.03 * h ** 1.72, 10) * log10(d_3d) -
+            _min(0.044 * h ** 1.72, 14.77) +
+            0.002 * log10(h) * d_3d
+            )
+    else:
+        PL1 = (
+                20 * log10(40 * M_PI * d_bp_3d * fc_ghz / 3) +
+                _min(0.03 * h ** 1.72, 10) * log10(d_bp_3d) -
+                _min(0.044 * h ** 1.72, 14.77) +
+                0.002 * log10(h) * d_bp_3d
+                )
+        PL_los = PL2 = PL1 + 40 * log10(d_3d / d_bp)
+
+    PL_nlos_prime = (
+        161.04 - 7.1 * log10(W) + 7.5 * log10(h) -
+        (24.37 - 3.7 * (h / h_bs) ** 2) * log10(h_bs) +
+        (43.42 - 3.1 * log10(h_bs)) * (log10(d_3d) - 3) +
+        20 * log10(fc_ghz) - (3.2 * (log10(11.75 * h_ue)) ** 2 - 4.97)
+        )
+
+    PL_nlos = _max(PL_los, PL_nlos_prime)
+
+    return (PL_los, PL_nlos)
+
+
+cdef (double, double) _urban_macro_losses(
+        double fc_ghz, double d_2d,
+        double h_bs, double h_ue,
+        double h_e,
+        ) nogil:
+    # all parameters (except freq) in meters;
+    # see 3GPP TR 38.901 version 16.1.0 Release 16;
+    # Table 7.4.1-1: Pathloss models; "UMa"
+    # Important: h_e parameter is random, but can be set to 1, if h_ue < 13 m
+    # (see Note 1 in Table 7.4.1-1)
+
+    cdef:
+        double d_3d, d_bp_prime
+        double PL1, PL2, PL_los, PL_nlos_prime, PL_nlos
+
+
+    d_3d = sqrt(d_2d ** 2 + (h_bs - h_ue) ** 2)
+    d_bp_prime = 4 *( h_bs - h_e) * (h_ue - h_e) * fc_ghz * 1e9 / C
+
+    if d_2d < 10 or d_2d > 5000:
+        return (NAN, NAN)
+
+    if d_2d <= d_bp_prime:
+        PL_los = PL1 = (
+            28. + 22 * log10(d_3d) + 20 * log10(fc_ghz)
+            )
+    else:
+        PL_los = PL2 = (
+            28. + 40 * log10(d_3d) + 20 * log10(fc_ghz) -
+            9 * log10(d_bp_prime ** 2 + (h_bs - h_ue) ** 2)
+            )
+
+    PL_nlos_prime = (
+        13.54 + 39.08 * log10(d_3d) + 20 * log10(fc_ghz) - 0.6 * (h_ue - 1.5)
+        )
+
+    PL_nlos = _max(PL_los, PL_nlos_prime)
+
+    return (PL_los, PL_nlos)
+
+
+cdef (double, double) _urban_micro_losses(
+        double fc_ghz, double d_2d,
+        double h_bs, double h_ue,
+        ) nogil:
+    # all parameters (except freq) in meters;
+    # see 3GPP TR 38.901 version 16.1.0 Release 16;
+    # Table 7.4.1-1: Pathloss models; "UMi - Street Canyon"
+
+    cdef:
+        double d_3d, d_bp_prime
+        double PL1, PL2, PL_los, PL_nlos_prime, PL_nlos
+        double h_e = 1.  # m
+
+
+    d_3d = sqrt(d_2d ** 2 + (h_bs - h_ue) ** 2)
+    d_bp_prime = 4 *( h_bs - h_e) * (h_ue - h_e) * fc_ghz * 1e9 / C
+
+    if d_2d < 10 or d_2d > 5000:
+        return (NAN, NAN)
+
+    if d_2d <= d_bp_prime:
+        PL_los = PL1 = (
+            32.4 + 21 * log10(d_3d) + 20 * log10(fc_ghz)
+            )
+    else:
+        PL_los = PL2 = (
+            32.4 + 40 * log10(d_3d) + 20 * log10(fc_ghz) -
+            9.5 * log10(d_bp_prime ** 2 + (h_bs - h_ue) ** 2)
+            )
+
+    PL_nlos_prime = (
+        22.4 + 35.3 * log10(d_3d) + 21.3 * log10(fc_ghz) - 0.3 * (h_ue - 1.5)
+        )
+
+    PL_nlos = _max(PL_los, PL_nlos_prime)
+
+    return (PL_los, PL_nlos)
+
+
+def rural_macro_losses_cython(
+        fc_ghz, d_2d,
+        h_bs, h_ue,
+        W, h,
+        out_PL_los=None,
+        out_PL_nlos=None
+        ):
+
+    cdef:
+
+        np.ndarray[double] _fc_ghz, _d_2d, _h_bs, _h_ue, _W, _h,
+        np.ndarray[double] _out_PL_los, _out_PL_nlos
+
+        int i, size
+
+    it = np.nditer(
+        [
+            fc_ghz, d_2d, h_bs, h_ue, W, h,
+            out_PL_los, out_PL_nlos
+            ],
+        flags=['external_loop', 'buffered', 'delay_bufalloc'],
+        op_flags=[
+            ['readonly'], ['readonly'], ['readonly'], ['readonly'],
+            ['readonly'], ['readonly'],
+            ['readwrite', 'allocate'], ['readwrite', 'allocate'],
+            ],
+        op_dtypes=[
+            'float64', 'float64', 'float64', 'float64', 'float64', 'float64',
+            'float64', 'float64'
+            ]
+        )
+
+    # it would be better to use the context manager but
+    # "with it:" requires numpy >= 1.14
+
+    it.reset()
+
+    for itup in it:
+        _fc_ghz = itup[0]
+        _d_2d = itup[1]
+        _h_bs = itup[2]
+        _h_ue = itup[3]
+        _W = itup[4]
+        _h = itup[5]
+        _out_PL_los = itup[6]
+        _out_PL_nlos = itup[7]
+
+        size = _fc_ghz.shape[0]
+
+        for i in prange(size, nogil=True):
+
+            (
+                _out_PL_los[i],
+                _out_PL_nlos[i],
+                ) = _rural_macro_losses(
+                    _fc_ghz[i],
+                    _d_2d[i],
+                    _h_bs[i],
+                    _h_ue[i],
+                    _W[i],
+                    _h[i],
+                    )
+
+    return it.operands[6:8]
+
+
+def urban_macro_losses_cython(
+        fc_ghz, d_2d,
+        h_bs, h_ue,
+        h_e,
+        out_PL_los=None,
+        out_PL_nlos=None
+        ):
+
+    cdef:
+
+        np.ndarray[double] _fc_ghz, _d_2d, _h_bs, _h_ue, _h_e,
+        np.ndarray[double] _out_PL_los, _out_PL_nlos
+
+        int i, size
+
+    it = np.nditer(
+        [
+            fc_ghz, d_2d, h_bs, h_ue, h_e,
+            out_PL_los, out_PL_nlos
+            ],
+        flags=['external_loop', 'buffered', 'delay_bufalloc'],
+        op_flags=[
+            ['readonly'], ['readonly'], ['readonly'], ['readonly'],
+            ['readonly'],
+            ['readwrite', 'allocate'], ['readwrite', 'allocate'],
+            ],
+        op_dtypes=[
+            'float64', 'float64', 'float64', 'float64', 'float64',
+            'float64', 'float64'
+            ]
+        )
+
+    # it would be better to use the context manager but
+    # "with it:" requires numpy >= 1.14
+
+    it.reset()
+
+    for itup in it:
+        _fc_ghz = itup[0]
+        _d_2d = itup[1]
+        _h_bs = itup[2]
+        _h_ue = itup[3]
+        _h_e = itup[4]
+        _out_PL_los = itup[5]
+        _out_PL_nlos = itup[6]
+
+        size = _fc_ghz.shape[0]
+
+        for i in prange(size, nogil=True):
+
+            (
+                _out_PL_los[i],
+                _out_PL_nlos[i],
+                ) = _urban_macro_losses(
+                    _fc_ghz[i],
+                    _d_2d[i],
+                    _h_bs[i],
+                    _h_ue[i],
+                    _h_e[i],
+                    )
+
+    return it.operands[5:7]
+
+
+def urban_micro_losses_cython(
+        fc_ghz, d_2d,
+        h_bs, h_ue,
+        out_PL_los=None,
+        out_PL_nlos=None
+        ):
+
+    cdef:
+
+        np.ndarray[double] _fc_ghz, _d_2d, _h_bs, _h_ue, _h_e,
+        np.ndarray[double] _out_PL_los, _out_PL_nlos
+
+        int i, size
+
+    it = np.nditer(
+        [
+            fc_ghz, d_2d, h_bs, h_ue,
+            out_PL_los, out_PL_nlos
+            ],
+        flags=['external_loop', 'buffered', 'delay_bufalloc'],
+        op_flags=[
+            ['readonly'], ['readonly'], ['readonly'], ['readonly'],
+            ['readwrite', 'allocate'], ['readwrite', 'allocate'],
+            ],
+        op_dtypes=[
+            'float64', 'float64', 'float64', 'float64',
+            'float64', 'float64'
+            ]
+        )
+
+    # it would be better to use the context manager but
+    # "with it:" requires numpy >= 1.14
+
+    it.reset()
+
+    for itup in it:
+        _fc_ghz = itup[0]
+        _d_2d = itup[1]
+        _h_bs = itup[2]
+        _h_ue = itup[3]
+        _out_PL_los = itup[4]
+        _out_PL_nlos = itup[5]
+
+        size = _fc_ghz.shape[0]
+
+        for i in prange(size, nogil=True):
+
+            (
+                _out_PL_los[i],
+                _out_PL_nlos[i],
+                ) = _urban_micro_losses(
+                    _fc_ghz[i],
+                    _d_2d[i],
+                    _h_bs[i],
+                    _h_ue[i],
+                    )
+
+    return it.operands[4:6]
diff --git a/pycraf/pathprof/helper.py b/pycraf/pathprof/helper.py
index cd02c72c0..bdcdc0182 100644
--- a/pycraf/pathprof/helper.py
+++ b/pycraf/pathprof/helper.py
@@ -485,7 +485,7 @@ def terrain_cmap_factory(sealevel=0.5, vmax=1200):
         >>> terrain_cmap, terrain_norm = terrain_cmap_factory(vmax=vmax)  # doctest: +SKIP
         >>> plt.imshow(  # doctest: +SKIP
         ...     heights, cmap=terrain_cmap, norm=terrain_norm,
-        ...     vmin=vmin, vmax=vmax
+        ...     # vmin=vmin, vmax=vmax  # deprecated in newer matplotlib versions
         ...     )
 
     Parameters
@@ -493,8 +493,9 @@ def terrain_cmap_factory(sealevel=0.5, vmax=1200):
     sealevel : float
         The sealevel value.
     vmax : float
-        Maximum height to cover in the colormap
-        (one should call plt.imshow with the same vmax option!)
+        Maximum height to cover in the colormap (Default: 1200)
+        (in older matplotlib versions, one should call plt.imshow with the
+        same `vmax` option, and `vmin=-20`!)
 
     Returns
     -------
diff --git a/pycraf/pathprof/imt.py b/pycraf/pathprof/imt.py
new file mode 100644
index 000000000..7722329fb
--- /dev/null
+++ b/pycraf/pathprof/imt.py
@@ -0,0 +1,624 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from __future__ import (
+    absolute_import, unicode_literals, division, print_function
+    )
+
+# from functools import partial, lru_cache
+from astropy import units as apu
+import numpy as np
+
+from . import cyimt
+from .. import conversions as cnv
+from .. import utils
+# import ipdb
+
+
+__all__ = [
+    'imt_rural_macro_losses', 'imt_urban_macro_losses',
+    'imt_urban_micro_losses', 'P_pusch', 'UE_MIN_MAX_DISTANCES',
+    'clutter_imt',
+    ]
+
+
+UE_MIN_MAX_DISTANCES = {
+    'rural_macro': (10, 10000),
+    'urban_macro': (10, 5000),
+    'urban_micro': (10, 5000),
+    }
+
+# UE_MIN_MAX_DISTANCES.__doc__ = '''
+#     Minimum and maximum distances that UEs have from base stations in each of
+#     the IMT propagation models. Below and beyond these, the path losses will
+#     be NaN.
+#     '''
+
+# Note, we have to curry the quantities here, because Cython produces
+# "built-in" functions that don't provide a signature (such that
+# ranged_quantity_input fails)
+
+
+def _Qinv(x):
+    # Note, this is *not* identical to cyprop._I_helper
+    # only good between 1.e-6 and 0.5
+    # See R-Rec P.1546
+
+    x = np.atleast_1d(x).copy()
+    mask = x > 0.5
+    x[mask] = 1 - x[mask]
+
+    T = np.sqrt(-2 * np.log(x))
+    Z = (
+        (
+            ((0.010328 * T + 0.802853) * T) + 2.515516698
+            ) /
+        (
+            ((0.001308 * T + 0.189269) * T + 1.432788) * T + 1.
+            )
+        )
+
+    Q = T - Z
+    Q[mask] *= -1
+    return Q
+
+
+# def Qinv(x):
+#     # larger x range than the approximation given in P.1546?
+#     # definitely much slower
+
+#     from scipy.stats import norm as qnorm
+
+#     x = np.atleast_1d(x).copy()
+
+#     mask = x > 0.5
+#     x[mask] = 1 - x[mask]
+
+#     Q = -qnorm.ppf(x, 0)
+#     Q[mask] *= -1
+
+#     return Q
+
+
+def _clutter_imt(
+        freq,
+        dist,
+        location_percent,
+        num_end_points=1,
+        ):
+
+    assert num_end_points in [1, 2]
+
+    L_l = 23.5 + 9.6 * np.log10(freq)
+    L_s = 32.98 + 23.9 * np.log10(dist) + 3.0 * np.log10(freq)
+
+    L_clutter = -5 * np.log10(
+        np.power(10, -0.2 * L_l) + np.power(10, -0.2 * L_s)
+        ) - 6 * _Qinv(location_percent / 100.)
+
+    if num_end_points == 2:
+        L_clutter *= 2
+
+    return L_clutter
+
+
+@utils.ranged_quantity_input(
+    freq=(2, 67, apu.GHz),
+    dist=(0.25, None, apu.km),
+    location_percent=(0, 100, apu.percent),
+    strip_input_units=True, output_unit=cnv.dB
+    )
+def clutter_imt(
+        freq,
+        dist,
+        location_percent,
+        num_end_points=1,
+        ):
+    '''
+    Calculate the Clutter loss according to IMT.CLUTTER document (method 2).
+
+    Parameters
+    ----------
+    freq : `~astropy.units.Quantity`
+        Frequency of radiation [GHz]
+    dist : `~astropy.units.Quantity`
+        Distance between Tx/Rx antennas [km]
+
+        Minimal distance must be 0.25 km (single endpoint clutter) or 1 km
+        (if both endpoints are to be corrected for clutter)
+    location_percent : `~astropy.units.Quantity`
+        Percentage of locations for which the clutter loss `L_clutter`
+        (calculated with this function) will not be exceeded [%]
+    num_end_points : int, optional
+        number of endpoints affected by clutter, allowed values: 1, 2
+
+    Returns
+    -------
+    L_clutter : `~astropy.units.Quantity`
+        Clutter loss [dB]
+
+    Notes
+    -----
+    - The algorithm is independent of effective antenna height (w.r.t.
+      clutter height), i.e., it doesn't distinguish between terminals which
+      are close to the ground and those closer to the top of the building.
+      However, the model is only appropriate if the terminal is "in the
+      clutter", below the rooftops.
+    - The result of this function is to be understood as a cumulative
+      value. For example, if `location_percent = 2%`, it means that for
+      2% of all possible locations, the clutter loss will not exceed the
+      returned `L_clutter` value, for the remaining 98% of locations it
+      will therefore be lower than `L_clutter`. The smaller `location_percent`,
+      the smaller the returned `L_clutter`, i.e., low clutter attenuations
+      are more unlikely.
+    - This model was proposed by ITU study group SG 3 to replace
+      `~pycraf.pathprof.clutter_correction` for IMT 5G studies (especially,
+      at higher frequencies, where multipath effects play a role in
+      urban and suburban areas).
+    '''
+
+    return _clutter_imt(
+        freq,
+        dist,
+        location_percent,
+        num_end_points=num_end_points,
+        )
+
+
+@utils.ranged_quantity_input(
+    freq=(0.5, 30, apu.GHz),
+    dist_2d=(0, 100000, apu.m),
+    h_bs=(10, 150, apu.m),
+    h_ue=(1, 10, apu.m),
+    W=(5, 50, apu.m),
+    h=(5, 50, apu.m),
+    strip_input_units=True, allow_none=False,
+    output_unit=(cnv.dB, cnv.dB, cnv.dimless)
+    )
+def imt_rural_macro_losses(
+        freq,
+        dist_2d,
+        h_bs=35 * apu.m, h_ue=1.5 * apu.m,
+        W=20 * apu.m, h=5 * apu.m,
+        ):
+    '''
+    Calculate los/non-los propagation losses Rural-Macro IMT scenario.
+
+    The computation is in accordance with `3GPP TR 38.901 Table 7.4.1-1
+    <https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=3173>`_
+
+    Parameters
+    ----------
+    freq : `~astropy.units.Quantity`
+        Frequency of radiation [GHz]
+    dist_2d : `~astropy.units.Quantity`
+        Distance on the ground between BS and UE device [m]
+        Note: Well-defined only for distances between 10 m and 10 km.
+    h_bs : `~astropy.units.Quantity`, optional
+        Basestation height [m] (Default: 35 m)
+    h_ue : `~astropy.units.Quantity`, optional
+        User equipment height [m] (Default: 1.5 m)
+    W : `~astropy.units.Quantity`, optional
+        Average street width [m] (Default: 20 m)
+    h : `~astropy.units.Quantity`, optional
+        Average building height [m] (Default: 5 m)
+
+    Returns
+    -------
+    PL_los : `~astropy.units.Quantity`
+        Path loss for line-of-sight cases [dB]
+    PL_nlos : `~astropy.units.Quantity`
+        Path loss for non-line-of-sight cases [dB]
+    los_prob : `~astropy.units.Quantity`
+        Probability for a path to be line-of-sight. [dimless]
+        (see Notes and Examples)
+
+    Notes
+    -----
+    - In statistical simulations, the LoS and Non-LoS cases occur with
+      certain probabilities. For sampling of path losses the return
+      parameter `los_prob` can be used, which accounts for the likelihoods
+      according to `3GPP TR 38.901 Table 7.4.2-1
+      <https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=3173>`_
+
+
+    Examples
+    --------
+
+    A typical usage, which also accounts for the line-of-sight probabilities,
+    would be::
+
+        >>> import numpy as np
+        >>> from pycraf import conversions as cnv
+        >>> from pycraf import pathprof
+        >>> from astropy import units as u
+        >>> from astropy.utils.misc import NumpyRNGContext
+
+        >>> freq = 1 * u.GHz
+        >>> h_bs, h_ue = 35 * u.m, 1.5 * u.m
+        >>> distances = [5, 20, 1000, 20000] * u.m  # 2D distances
+        >>> # Note: too small or large distances will lead to NaN values
+
+        >>> PL_los, PL_nlos, los_prob = pathprof.imt_rural_macro_losses(
+        ...     freq, distances, h_bs=h_bs, h_ue=h_ue
+        ...     )
+        >>> PL_los  # doctest: +FLOAT_CMP
+        <Decibel [        nan, 64.38071083, 94.57828524,         nan] dB>
+        >>> PL_nlos  # doctest: +FLOAT_CMP
+        <Decibel [         nan,  65.10847815, 119.54294519,          nan] dB>
+        >>> los_prob  # doctest: +FLOAT_CMP
+        <Quantity [1.00000000e+00, 9.90049834e-01, 3.71576691e-01, 2.08186856e-09]>
+
+        >>> # randomly assign LOS or Non-LOS type to UE (according to above prob.)
+        >>> with NumpyRNGContext(0):
+        ...     los_type = np.random.uniform(0, 1, distances.size) < los_prob
+        >>> # note: los_type == True : LOS
+        >>> #       los_type == False: Non-LOS
+        >>> los_type
+        array([ True,  True, False, False])
+        >>> PL = np.where(
+        ...     los_type, PL_los.to_value(cnv.dB), PL_nlos.to_value(cnv.dB)
+        ...     ) * cnv.dB
+        >>> PL  # doctest: +FLOAT_CMP
+        <Decibel [         nan,  64.38071083, 119.54294519,          nan] dB>
+
+    '''
+
+    pl_los, pl_nlos = cyimt.rural_macro_losses_cython(
+        freq, dist_2d, h_bs, h_ue, W, h,
+        )
+
+    # probability to have a LOS path;
+    # see 3GPP TR 38.901, Table 7.4.2
+    los_prob = np.exp(-(dist_2d - 10) / 1000)
+    los_prob[dist_2d < 10] = 1.
+    los_prob = np.broadcast_to(los_prob, pl_los.shape)
+
+    return pl_los, pl_nlos, los_prob
+
+
+@utils.ranged_quantity_input(
+    freq=(0.5, 30, apu.GHz),
+    dist_2d=(0, 100000, apu.m),
+    h_bs=(10, 150, apu.m),  # according to 3GPP TR 38.901 only 25 m allowed!?
+    h_ue=(1, 13, apu.m),  # for larger h_ue, need to implement "C" function
+    h_e=(0, 20, apu.m),
+    strip_input_units=True, allow_none=False,
+    output_unit=(cnv.dB, cnv.dB, cnv.dimless)
+    )
+def imt_urban_macro_losses(
+        freq,
+        dist_2d,
+        h_bs=25 * apu.m, h_ue=1.5 * apu.m,
+        # h_e=1 * apu.m,
+        ):
+    '''
+    Calculate los/non-los propagation losses Rural-Macro IMT scenario.
+
+    The computation is in accordance with `3GPP TR 38.901 Table 7.4.1-1
+    <https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=3173>`_
+
+    Parameters
+    ----------
+    freq : `~astropy.units.Quantity`
+        Frequency of radiation [GHz]
+    dist_2d : `~astropy.units.Quantity`
+        Distance on the ground between BS and UE device [m]
+        Note: Well-defined only for distances between 10 m and 5 km.
+    h_bs : `~astropy.units.Quantity`, optional
+        Basestation height [m] (Default: 35 m)
+    h_ue : `~astropy.units.Quantity`, optional
+        User equipment height [m] (Default: 1.5 m)
+        Note: in the `pycraf` implementation this is restricted to
+        `h_ue < 13 m`. 3GPP TR 38.901 also has a model for heights up
+        to 22.5 m, but this involves a `h_e` different from 1 m and is
+        probabilistic, which would make the interface much more complicated.
+
+    Returns
+    -------
+    PL_los : `~astropy.units.Quantity`
+        Path loss for line-of-sight cases [dB]
+    PL_nlos : `~astropy.units.Quantity`
+        Path loss for non-line-of-sight cases [dB]
+    los_prob : `~astropy.units.Quantity`
+        Probability for a path to be line-of-sight. [dimless]
+        (see Notes and Examples)
+
+    Notes
+    -----
+    - In statistical simulations, the LoS and Non-LoS cases occur with
+      certain probabilities. For sampling of path losses the return
+      parameter `los_prob` can be used, which accounts for the likelihoods
+      according to `3GPP TR 38.901 Table 7.4.2-1
+      <https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=3173>`_
+
+
+    Examples
+    --------
+
+    A typical usage, which also accounts for the line-of-sight probabilities,
+    would be::
+
+        >>> import numpy as np
+        >>> from pycraf import conversions as cnv
+        >>> from pycraf import pathprof
+        >>> from astropy import units as u
+        >>> from astropy.utils.misc import NumpyRNGContext
+
+        >>> freq = 1 * u.GHz
+        >>> h_bs, h_ue = 25 * u.m, 1.5 * u.m
+        >>> distances = [5, 20, 1000, 20000] * u.m  # 2D distances
+        >>> # Note: too small or large distances will lead to NaN values
+
+        >>> PL_los, PL_nlos, los_prob = pathprof.imt_urban_macro_losses(
+        ...     freq, distances, h_bs=h_bs, h_ue=h_ue
+        ...     )
+        >>> PL_los  # doctest: +FLOAT_CMP
+        <Decibel [         nan,  60.76626079, 108.24721393,          nan] dB>
+        >>> PL_nlos  # doctest: +FLOAT_CMP
+        <Decibel [         nan,  71.74479418, 130.78468516,          nan] dB>
+        >>> los_prob  # doctest: +FLOAT_CMP
+        <Quantity [1.00000000e+00, 9.72799557e-01, 1.80001255e-02, 9.00000000e-04]>
+
+        >>> # randomly assign LOS or Non-LOS type to UE (according to above prob.)
+        >>> with NumpyRNGContext(0):
+        ...     los_type = np.random.uniform(0, 1, distances.size) < los_prob
+        >>> # note: los_type == True : LOS
+        >>> #       los_type == False: Non-LOS
+        >>> los_type
+        array([ True,  True, False, False])
+        >>> PL = np.where(
+        ...     los_type, PL_los.to_value(cnv.dB), PL_nlos.to_value(cnv.dB)
+        ...     ) * cnv.dB
+        >>> PL  # doctest: +FLOAT_CMP
+        <Decibel [         nan,  60.76626079, 130.78468516,          nan] dB>
+
+    '''
+    # for h_ue > 13 m, need to implement "C" function; then "h_e" is needed!
+    # h_e : `~astropy.units.Quantity`, optional
+    #     Effective environment height [m] (Default: 1 m)
+    #     Important: for `h_ue > 13 m`, this is subject to random sampling;
+    #     see `3GPP TR 38.901 Table 7.4.1-1 Note 1
+    #     <https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=3173>`_
+
+    h_e = 1.
+    pl_los, pl_nlos = cyimt.urban_macro_losses_cython(
+        freq, dist_2d, h_bs, h_ue, h_e,
+        )
+
+    # probability to have a LOS path;
+    # see 3GPP TR 38.901, Table 7.4.2
+    _x = 18 / dist_2d
+    los_prob = (
+        # for h_e > 13 m, a correction factor would be necessary!!!
+        _x + np.exp(-18 / 63 / _x) * (1 - _x)
+        )
+    los_prob[dist_2d < 18] = 1.
+    los_prob = np.broadcast_to(los_prob, pl_los.shape)
+
+    return pl_los, pl_nlos, los_prob
+
+
+@utils.ranged_quantity_input(
+    freq=(0.5, 100, apu.GHz),
+    dist_2d=(0, 100000, apu.m),
+    h_bs=(10, 150, apu.m),  # according to 3GPP TR 38.901 only 10 m allowed!?
+    h_ue=(1, 22.5, apu.m),
+    strip_input_units=True, allow_none=False,
+    output_unit=(cnv.dB, cnv.dB, cnv.dimless)
+    )
+def imt_urban_micro_losses(
+        freq,
+        dist_2d,
+        h_bs=10 * apu.m, h_ue=1.5 * apu.m,
+        ):
+    '''
+    Calculate los/non-los propagation losses Rural-Macro IMT scenario.
+
+    The computation is in accordance with `3GPP TR 38.901 Table 7.4.1-1
+    <https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=3173>`_
+
+    Parameters
+    ----------
+    freq : `~astropy.units.Quantity`
+        Frequency of radiation [GHz]
+    dist_2d : `~astropy.units.Quantity`
+        Distance on the ground between BS and UE device [m]
+        Note: Well-defined only for distances between 10 m and 5 km.
+    h_bs : `~astropy.units.Quantity`, optional
+        Basestation height [m] (Default: 35 m)
+    h_ue : `~astropy.units.Quantity`, optional
+        User equipment height [m] (Default: 1.5 m)
+
+    Returns
+    -------
+    PL_los : `~astropy.units.Quantity`
+        Path loss for line-of-sight cases [dB]
+    PL_nlos : `~astropy.units.Quantity`
+        Path loss for non-line-of-sight cases [dB]
+    los_prob : `~astropy.units.Quantity`
+        Probability for a path to be line-of-sight. [dimless]
+        (see Notes and Examples)
+
+    Notes
+    -----
+    - In statistical simulations, the LoS and Non-LoS cases occur with
+      certain probabilities. For sampling of path losses the return
+      parameter `los_prob` can be used, which accounts for the likelihoods
+      according to `3GPP TR 38.901 Table 7.4.2-1
+      <https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=3173>`_
+
+
+    Examples
+    --------
+
+    A typical usage, which also accounts for the line-of-sight probabilities,
+    would be::
+
+        >>> import numpy as np
+        >>> from pycraf import conversions as cnv
+        >>> from pycraf import pathprof
+        >>> from astropy import units as u
+        >>> from astropy.utils.misc import NumpyRNGContext
+
+        >>> freq = 1 * u.GHz
+        >>> h_bs, h_ue = 10 * u.m, 1.5 * u.m
+        >>> distances = [5, 20, 1000, 20000] * u.m  # 2D distances
+        >>> # Note: too small or large distances will lead to NaN values
+
+        >>> PL_los, PL_nlos, los_prob = pathprof.imt_urban_micro_losses(
+        ...     freq, distances, h_bs=h_bs, h_ue=h_ue
+        ...     )
+        >>> PL_los  # doctest: +FLOAT_CMP
+        <Decibel [         nan,  60.47880565, 118.53377126,          nan] dB>
+        >>> PL_nlos  # doctest: +FLOAT_CMP
+        <Decibel [         nan,  69.59913521, 128.3005538 ,          nan] dB>
+        >>> los_prob  # doctest: +FLOAT_CMP
+        <Quantity [1.00000000e+00, 9.57375342e-01, 1.80000000e-02, 9.00000000e-04]>
+
+        >>> # randomly assign LOS or Non-LOS type to UE (according to above prob.)
+        >>> with NumpyRNGContext(0):
+        ...     los_type = np.random.uniform(0, 1, distances.size) < los_prob
+        >>> # note: los_type == True : LOS
+        >>> #       los_type == False: Non-LOS
+        >>> los_type
+        array([ True,  True, False, False])
+        >>> PL = np.where(
+        ...     los_type, PL_los.to_value(cnv.dB), PL_nlos.to_value(cnv.dB)
+        ...     ) * cnv.dB
+        >>> PL  # doctest: +FLOAT_CMP
+        <Decibel [         nan,  60.47880565, 128.3005538 ,          nan] dB>
+
+    '''
+
+    pl_los, pl_nlos = cyimt.urban_micro_losses_cython(
+        freq, dist_2d, h_bs, h_ue,
+        )
+
+    # probability to have a LOS path;
+    # see 3GPP TR 38.901, Table 7.4.2
+    _x = 18 / dist_2d
+    los_prob = (
+        _x + np.exp(-0.5 / _x) * (1 - _x)
+        )
+    los_prob[dist_2d < 18] = 1.
+    los_prob = np.broadcast_to(los_prob, pl_los.shape)
+
+    return pl_los, pl_nlos, los_prob
+
+
+@utils.ranged_quantity_input(
+    P_cmax=(-20, 50, cnv.dBm),  # range pretty large, but what is realistic?
+    P_0_pusch=(-200, 100, cnv.dBm),  # dito
+    alpha=(0, 1, cnv.dimless),
+    PL=(-100, 300, cnv.dB),
+    strip_input_units=True, allow_none=False,
+    output_unit=(cnv.dBm)
+    )
+def P_pusch(P_cmax, M_pusch, P_0_pusch, alpha, PL):
+    '''
+    Calculate power output level after UE power control.
+
+    See `ITU-R Rec. M.2101-0 <https://www.itu.int/rec/R-REC-M.2101/en>`_
+    Section 4.1.
+
+    Parameters
+    ----------
+    P_cmax : `~astropy.units.Quantity`
+        Maximum transmit power [dBm]
+    M_pusch : `numpy.ndarray`, int
+        Number of allocated resource blocks (RBs)
+
+        Is this the bandwidth per carrier divided by the RB bandwidth (
+        typically 180 kHz) and number of UE devices associated to each
+        carrier?
+    P_0_pusch : `numpy.ndarray`
+        Initial receive target UE power level per RB [dBm]
+    alpha : `~astropy.units.Quantity`
+        Balancing factor for UEs with bad channel and UEs with good channel
+    PL : `~astropy.units.Quantity`
+        Path loss between UE and its associated BS [dB]
+        One should use one of the functions
+        `~pathprof.imt_rural_macro_losses`,
+        `~pathprof.imt_urban_macro_losses`, or
+        `~pathprof.imt_urban_micro_losses` to determine PL for the required
+        scenario.
+
+        Note: Antenna gains should be included, so this is rather the
+        coupling loss than the path loss, unlike what is stated in
+        `ITU-R Rec. M.2101-0 <https://www.itu.int/rec/R-REC-M.2101/en>`_.
+        ``
+
+    Returns
+    -------
+    P_pusch : `~astropy.units.Quantity`
+        Transmit power of the terminal [dBm]
+
+
+
+    Examples
+    --------
+
+    A typical usage would be::
+
+        >>> import numpy as np
+        >>> from pycraf import conversions as cnv
+        >>> from pycraf import pathprof
+        >>> from astropy import units as u
+        >>> from astropy.utils.misc import NumpyRNGContext
+
+        >>> freq = 6.65 * u.GHz
+        >>> h_bs, h_ue = 10 * u.m, 1.5 * u.m
+        >>> distances = [20, 100, 500, 1000] * u.m  # 2D distances
+
+        >>> PL_los, PL_nlos, los_prob = pathprof.imt_urban_micro_losses(
+        ...     freq, distances, h_bs=h_bs, h_ue=h_ue
+        ...     )
+
+        >>> # randomly assign LOS or Non-LOS type to UE (according to above prob.)
+        >>> with NumpyRNGContext(0):
+        ...     los_type = np.random.uniform(0, 1, distances.size) < los_prob
+        >>> # note: los_type == True : LOS
+        >>> #       los_type == False: Non-LOS
+        >>> PL = np.where(
+        ...     los_type, PL_los.to_value(cnv.dB), PL_nlos.to_value(cnv.dB)
+        ...     ) * cnv.dB
+        >>> PL  # doctest: +FLOAT_CMP
+        <Decibel [ 76.93523856, 110.58128371, 135.20195715, 145.82665484] dB>
+
+        >>> # Assume some antenna gains:
+        >>> G_bs, G_ue = 20 * cnv.dBi, 5 * cnv.dBi
+        >>> CL = (
+        ...     PL.to_value(cnv.dB) -
+        ...     G_bs.to_value(cnv.dBi) -
+        ...     G_ue.to_value(cnv.dBi)
+        ...     ) * cnv.dB
+
+
+        >>> bw_carrier = 100 * u.MHz
+        >>> bw_rb = 180 * u.kHz  # resource block bandwidth
+        >>> num_ue = 3  # 3 UEs per BS sector and carrier
+        >>> M_pusch = int(bw_carrier / bw_rb / num_ue)
+        >>> M_pusch
+        185
+
+        >>> P_cmax = 23 * cnv.dBm
+        >>> P_0_pusch = -95.5 * cnv.dBm
+        >>> alpha = 0.8 * cnv.dimless
+
+        >>> P_pusch = pathprof.P_pusch(
+        ...     P_cmax, M_pusch, P_0_pusch, alpha, CL
+        ...     )
+        >>> P_pusch.to(cnv.dBm)
+        <Decibel [-31.28009187,  -4.36325575,  15.333283  ,  23.          ] dB(mW)>
+
+    '''
+
+    P_pusch = 10 * np.log10(M_pusch) + P_0_pusch + alpha * PL
+    P_pusch[P_pusch > P_cmax] = P_cmax
+
+    return P_pusch
+
+
+if __name__ == '__main__':
+    print('This not a standalone python program! Use as module.')
diff --git a/pycraf/pathprof/propagation.py b/pycraf/pathprof/propagation.py
index 115b01824..a3f20155d 100644
--- a/pycraf/pathprof/propagation.py
+++ b/pycraf/pathprof/propagation.py
@@ -22,7 +22,7 @@
     'PathProp',
     'loss_freespace', 'loss_troposcatter', 'loss_ducting',
     'loss_diffraction', 'loss_complete',
-    'clutter_correction', 'clutter_imt',
+    'clutter_correction',
     'height_map_data', 'atten_map_fast',
     'height_path_data', 'height_path_data_generic', 'atten_path_fast',
     'losses_complete',
@@ -472,132 +472,6 @@ def clutter_correction(
     return cyprop.clutter_correction_cython(h_g, zone, freq)
 
 
-def _Qinv(x):
-    # Note, this is *not* identical to cyprop._I_helper
-    # only good between 1.e-6 and 0.5
-    # See R-Rec P.1546
-
-    x = np.atleast_1d(x).copy()
-    mask = x > 0.5
-    x[mask] = 1 - x[mask]
-
-    T = np.sqrt(-2 * np.log(x))
-    Z = (
-        (
-            ((0.010328 * T + 0.802853) * T) + 2.515516698
-            ) /
-        (
-            ((0.001308 * T + 0.189269) * T + 1.432788) * T + 1.
-            )
-        )
-
-    Q = T - Z
-    Q[mask] *= -1
-    return Q
-
-
-# def Qinv(x):
-#     # larger x range than the approximation given in P.1546?
-#     # definitely much slower
-
-#     from scipy.stats import norm as qnorm
-
-#     x = np.atleast_1d(x).copy()
-
-#     mask = x > 0.5
-#     x[mask] = 1 - x[mask]
-
-#     Q = -qnorm.ppf(x, 0)
-#     Q[mask] *= -1
-
-#     return Q
-
-
-def _clutter_imt(
-        freq,
-        dist,
-        location_percent,
-        num_end_points=1,
-        ):
-
-    assert num_end_points in [1, 2]
-
-    L_l = 23.5 + 9.6 * np.log10(freq)
-    L_s = 32.98 + 23.9 * np.log10(dist) + 3.0 * np.log10(freq)
-
-    L_clutter = -5 * np.log10(
-        np.power(10, -0.2 * L_l) + np.power(10, -0.2 * L_s)
-        ) - 6 * _Qinv(location_percent / 100.)
-
-    if num_end_points == 2:
-        L_clutter *= 2
-
-    return L_clutter
-
-
-@utils.ranged_quantity_input(
-    freq=(2, 67, apu.GHz),
-    dist=(0.25, None, apu.km),
-    location_percent=(0, 100, apu.percent),
-    strip_input_units=True, output_unit=cnv.dB
-    )
-def clutter_imt(
-        freq,
-        dist,
-        location_percent,
-        num_end_points=1,
-        ):
-    '''
-    Calculate the Clutter loss according to IMT.CLUTTER document (method 2).
-
-    Parameters
-    ----------
-    freq : `~astropy.units.Quantity`
-        Frequency of radiation [GHz]
-    dist : `~astropy.units.Quantity`
-        Distance between Tx/Rx antennas [km]
-
-        Minimal distance must be 0.25 km (single endpoint clutter) or 1 km
-        (if both endpoints are to be corrected for clutter)
-    location_percent : `~astropy.units.Quantity`
-        Percentage of locations for which the clutter loss `L_clutter`
-        (calculated with this function) will not be exceeded [%]
-    num_end_points : int, optional
-        number of endpoints affected by clutter, allowed values: 1, 2
-
-    Returns
-    -------
-    L_clutter : `~astropy.units.Quantity`
-        Clutter loss [dB]
-
-    Notes
-    -----
-    - The algorithm is independent of effective antenna height (w.r.t.
-      clutter height), i.e., it doesn't distinguish between terminals which
-      are close to the ground and those closer to the top of the building.
-      However, the model is only appropriate if the terminal is "in the
-      clutter", below the rooftops.
-    - The result of this function is to be understood as a cumulative
-      value. For example, if `location_percent = 2%`, it means that for
-      2% of all possible locations, the clutter loss will not exceed the
-      returned `L_clutter` value, for the remaining 98% of locations it
-      will therefore be lower than `L_clutter`. The smaller `location_percent`,
-      the smaller the returned `L_clutter`, i.e., low clutter attenuations
-      are more unlikely.
-    - This model was proposed by ITU study group SG 3 to replace
-      `~pycraf.pathprof.clutter_correction` for IMT 5G studies (especially,
-      at higher frequencies, where multipath effects play a role in
-      urban and suburban areas).
-    '''
-
-    return _clutter_imt(
-        freq,
-        dist,
-        location_percent,
-        num_end_points=num_end_points,
-        )
-
-
 # TODO: do we want to convert output dictionary arrays to quantities?
 @utils.ranged_quantity_input(
     lon_t=(-180, 180, apu.deg),
diff --git a/pycraf/pathprof/setup_package.py b/pycraf/pathprof/setup_package.py
index 08d9ca1db..b419a6650 100644
--- a/pycraf/pathprof/setup_package.py
+++ b/pycraf/pathprof/setup_package.py
@@ -51,6 +51,12 @@ def get_extensions():
         comp_args['extra_compile_args'] = extra_compile_args
         # comp_args['extra_link_args'].append('-lgomp')
 
+    ext_module_pathprof_cyimt = Extension(
+        name='pycraf.pathprof.cyimt',
+        sources=[os.path.join(PYXDIR, 'cyimt.pyx')],
+        **comp_args
+        )
+
     ext_module_pathprof_cyprop = Extension(
         name='pycraf.pathprof.cyprop',
         sources=[os.path.join(PYXDIR, 'cyprop.pyx')],
@@ -63,4 +69,8 @@ def get_extensions():
         **comp_args
         )
 
-    return [ext_module_pathprof_cyprop, ext_module_pathprof_geodesics]
+    return [
+        ext_module_pathprof_cyimt,
+        ext_module_pathprof_cyprop,
+        ext_module_pathprof_geodesics
+        ]
diff --git a/setup.cfg b/setup.cfg
index 485c8e213..33f3ae6cf 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -56,7 +56,7 @@ github_project = bwinkel/pycraf
 # install_requires = astropy, scipy, matplotlib
 install_requires = astropy, scipy, matplotlib
 # version should be PEP440 compatible (https://www.python.org/dev/peps/pep-0440/)
-version = 1.1.0
+version = 1.1.1dev
 minimum_python_version = 3.6
 
 [entry_points]


Test patch:
diff --git a/pycraf/antenna/tests/test_antenna.py b/pycraf/antenna/tests/test_antenna.py
index 8fbe7b3c4..d09f165e9 100644
--- a/pycraf/antenna/tests/test_antenna.py
+++ b/pycraf/antenna/tests/test_antenna.py
@@ -460,6 +460,102 @@ def test_imt_advanced_sectoral_peak_sidelobe_pattern_400_to_6000_mhz():
         )
 
 
+
+
+def test_imt_advanced_sectoral_avg_sidelobe_pattern_400_to_6000_mhz():
+
+    _gfunc = imt.imt_advanced_sectoral_avg_sidelobe_pattern_400_to_6000_mhz
+    args_list = [
+        (-180, 180, apu.deg),
+        (-90, 90, apu.deg),
+        (None, None, cnv.dB),
+        (0, None, apu.deg),
+        (0, None, apu.deg),
+        (0, None, cnv.dimless),
+        (0, None, cnv.dimless),
+        (0, None, cnv.dimless),
+        (-90, 90, apu.deg),
+        (-90, 90, apu.deg),
+        ]
+    check_astro_quantities(_gfunc, args_list)
+
+    azims = np.arange(-180, 180, 0.5)[50:-50:150] * apu.deg
+    elevs = np.arange(-90, 90, 0.5)[30:-30:70] * apu.deg
+
+    G0 = 18. * cnv.dB
+    phi_3db = 65. * apu.deg
+    # theta_3db can be inferred in the following way:
+    theta_3db = 31000 / G0.to(cnv.dimless) / phi_3db.value * apu.deg
+    k_a, k_h, k_v = (0.7, 0.7, 0.3) * cnv.dimless
+    tilt_m, tilt_e = (0, 0) * apu.deg
+
+    bs_gain = _gfunc(
+        0 * apu.deg, 0 * apu.deg,
+        G0, phi_3db, theta_3db,
+        k_a, k_h, k_v,
+        tilt_m=tilt_m, tilt_e=tilt_e,
+        )
+
+    assert_quantity_allclose(bs_gain, G0)
+
+    bs_gains = _gfunc(
+        azims[np.newaxis], elevs[:, np.newaxis],
+        G0, phi_3db, theta_3db,
+        k_a, k_h, k_v,
+        tilt_m=tilt_m, tilt_e=tilt_e,
+        )
+    assert_quantity_allclose(
+        bs_gains,
+        [
+            [-9.45692316, -8.73264105, -7.99973768, -8.59935781, -9.45692316],
+            [-9.45692316, -6.23545959, -2.9756504 , -5.64264212, -9.45692316],
+            [-9.45692316,  1.55191039, 12.6917839 ,  3.57776872, -9.45692316],
+            [-9.45692316, -5.11617053, -0.72374921, -4.3173802 , -9.45692316],
+            [-9.45692316, -8.16416531, -6.85601954, -7.9262705 , -9.45692316],
+            ] * cnv.dB
+        )
+
+    tilt_m, tilt_e = (40, 0) * apu.deg
+
+    bs_gains = _gfunc(
+        azims[np.newaxis], elevs[:, np.newaxis],
+        G0, phi_3db, theta_3db,
+        k_a, k_h, k_v,
+        tilt_m=tilt_m, tilt_e=tilt_e,
+        )
+
+    assert_quantity_allclose(
+        bs_gains,
+        [
+            [-6.65167057, -4.21643061, -1.90092642, -3.78461285, -6.46818041],
+            [-9.04857002, -2.88936373, 17.95595854, -1.93430988, -8.57901892],
+            [-9.45692316,  3.69272869, -1.87877821, -2.04640382, -9.45692316],
+            [-9.45692316, -7.09347806, -7.4361134 , -6.90931547, -9.45692316],
+            [-9.45692316, -9.45692316, -9.45692316, -9.45692316, -9.45692316],
+            ] * cnv.dB
+        )
+
+    tilt_m, tilt_e = (0, 40) * apu.deg
+
+    bs_gains = _gfunc(
+        azims[np.newaxis], elevs[:, np.newaxis],
+        G0, phi_3db, theta_3db,
+        k_a, k_h, k_v,
+        tilt_m=tilt_m, tilt_e=tilt_e,
+        )
+
+    assert_quantity_allclose(
+        bs_gains,
+        [
+            [-9.45692316, -8.04001291, -6.60623694, -7.77927144, -9.45692316],
+            [-9.45692316,  4.15502284, 17.92899408,  6.65990894, -9.45692316],
+            [-9.45692316, -4.88763097, -0.26394973, -4.04678451, -9.45692316],
+            [-9.45692316, -6.99775686, -4.5093187 , -6.54521822, -9.45692316],
+            [-9.45692316, -8.60848865, -7.74995508, -8.45235874, -9.45692316],
+            ] * cnv.dB
+        )
+
+
 def test_fl_pattern():
 
     args_list = [
diff --git a/pycraf/pathprof/tests/test_propagation.py b/pycraf/pathprof/tests/test_propagation.py
index 80a006c10..c154d09da 100644
--- a/pycraf/pathprof/tests/test_propagation.py
+++ b/pycraf/pathprof/tests/test_propagation.py
@@ -1143,46 +1143,3 @@ def test_clutter_correction():
             loss,
             )
 
-
-def test_clutter_imt():
-
-    args_list = [
-        (2, 67, apu.GHz),
-        (0.25, None, apu.km),
-        (0, 100, apu.percent),
-        ]
-
-    check_astro_quantities(pathprof.clutter_imt, args_list)
-
-    freq = [2.5, 10, 20] * apu.GHz
-    dist = np.logspace(2.5, 4, 4)[:, np.newaxis] * apu.m
-
-    assert_quantity_allclose(
-        pathprof.clutter_imt(freq, dist, 2 * apu.percent, num_end_points=1),
-        [
-            [9.70031407, 11.67179752, 12.59465169],
-            [14.90452998, 20.26362512, 22.56060172],
-            [14.99471557, 20.77252115, 23.65890113],
-            [14.99509074, 20.77485862, 23.66473226]
-            ] * cnv.dB,
-        )
-
-    assert_quantity_allclose(
-        pathprof.clutter_imt(freq, dist, 50 * apu.percent, num_end_points=2),
-        [
-            [44.05089174, 47.99385864, 49.83956699],
-            [54.45932356, 65.17751385, 69.77146704],
-            [54.63969475, 66.19530591, 71.96806586],
-            [54.64044508, 66.19998085, 71.97972813]
-            ] * cnv.dB,
-        )
-
-    assert_quantity_allclose(
-        pathprof.clutter_imt(freq, dist, 90 * apu.percent, num_end_points=1),
-        [
-            [29.71581878, 31.68730223, 32.61015641],
-            [34.92003469, 40.27912984, 42.57610643],
-            [35.01022029, 40.78802587, 43.67440584],
-            [35.01059545, 40.79036334, 43.68023698]
-            ] * cnv.dB,
-        )
diff --git a/pycraf/pathprof/tests/test_propagation_imt.py b/pycraf/pathprof/tests/test_propagation_imt.py
new file mode 100644
index 000000000..bbbe5f370
--- /dev/null
+++ b/pycraf/pathprof/tests/test_propagation_imt.py
@@ -0,0 +1,252 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# from __future__ import absolute_import
+# from __future__ import division
+# from __future__ import print_function
+# from __future__ import unicode_literals
+
+import os
+import pytest
+# from functools import partial
+import numpy as np
+# from zipfile import ZipFile
+from numpy.testing import assert_equal, assert_allclose
+from astropy.tests.helper import assert_quantity_allclose, remote_data
+from astropy import units as apu
+from astropy.units import Quantity
+from ... import conversions as cnv
+from ... import pathprof
+from ...utils import check_astro_quantities
+# from astropy.utils.data import get_pkg_data_filename
+from astropy.utils.misc import NumpyRNGContext
+# import json
+# from itertools import product
+# import importlib
+
+
+TOL_KWARGS = {'atol': 1.e-4, 'rtol': 1.e-4}
+
+
+def test_clutter_imt():
+
+    args_list = [
+        (2, 67, apu.GHz),
+        (0.25, None, apu.km),
+        (0, 100, apu.percent),
+        ]
+
+    check_astro_quantities(pathprof.clutter_imt, args_list)
+
+    freq = [2.5, 10, 20] * apu.GHz
+    dist = np.logspace(2.5, 4, 4)[:, np.newaxis] * apu.m
+
+    assert_quantity_allclose(
+        pathprof.clutter_imt(freq, dist, 2 * apu.percent, num_end_points=1),
+        [
+            [9.70031407, 11.67179752, 12.59465169],
+            [14.90452998, 20.26362512, 22.56060172],
+            [14.99471557, 20.77252115, 23.65890113],
+            [14.99509074, 20.77485862, 23.66473226]
+            ] * cnv.dB,
+        )
+
+    assert_quantity_allclose(
+        pathprof.clutter_imt(freq, dist, 50 * apu.percent, num_end_points=2),
+        [
+            [44.05089174, 47.99385864, 49.83956699],
+            [54.45932356, 65.17751385, 69.77146704],
+            [54.63969475, 66.19530591, 71.96806586],
+            [54.64044508, 66.19998085, 71.97972813]
+            ] * cnv.dB,
+        )
+
+    assert_quantity_allclose(
+        pathprof.clutter_imt(freq, dist, 90 * apu.percent, num_end_points=1),
+        [
+            [29.71581878, 31.68730223, 32.61015641],
+            [34.92003469, 40.27912984, 42.57610643],
+            [35.01022029, 40.78802587, 43.67440584],
+            [35.01059545, 40.79036334, 43.68023698]
+            ] * cnv.dB,
+        )
+
+
+def test_imt_rural_macro_losses():
+
+    _func = pathprof.imt_rural_macro_losses
+    args_list = [
+        (0.5, 30, apu.GHz),
+        (0, 100000, apu.m),
+        (10, 150, apu.m),
+        (1, 10, apu.m),
+        (5, 50, apu.m),
+        (5, 50, apu.m),
+        ]
+
+    check_astro_quantities(_func, args_list)
+
+    freq = [1, 5] * apu.GHz
+    dist = np.array([5, 20, 1000, 20000])[:, np.newaxis] * apu.m
+
+    PL_los, PL_nlos, los_prob = _func(freq, dist)
+    print(PL_los)
+    print(PL_nlos)
+    print(los_prob)
+    assert_quantity_allclose(
+        PL_los,
+        [
+            [      np.nan,       np.nan],
+            [ 64.38071083,  78.36011092],
+            [ 94.57828524, 108.55768533],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )
+    assert_quantity_allclose(
+        PL_nlos,
+        [
+            [      np.nan,       np.nan],
+            [ 65.10847815,  79.08787824],
+            [119.54294519, 133.52234528],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )
+    assert_quantity_allclose(
+        los_prob,
+        [
+            [1.00000000e+00, 1.00000000e+00],
+            [9.90049834e-01, 9.90049834e-01],
+            [3.71576691e-01, 3.71576691e-01],
+            [2.08186856e-09, 2.08186856e-09],
+            ] * cnv.dimless,
+        )
+
+    PL_los, PL_nlos, los_prob = _func(freq, dist, h_bs=20 * apu.m)
+    print(PL_los)
+    print(PL_nlos)
+    assert_quantity_allclose(
+        PL_los,
+        [
+            [      np.nan,       np.nan],
+            [ 61.17035763,  75.14975772],
+            [ 98.00008131, 108.55367443],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )
+    assert_quantity_allclose(
+        PL_nlos,
+        [
+            [      np.nan,       np.nan],
+            [ 64.01112576,  77.99052585],
+            [125.64356998, 139.62297006],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )
+
+
+def test_imt_urban_macro_losses():
+
+    _func = pathprof.imt_urban_macro_losses
+    args_list = [
+        (0.5, 30, apu.GHz),
+        (0, 100000, apu.m),
+        (10, 150, apu.m),
+        (1, 13, apu.m),
+        ]
+
+    check_astro_quantities(_func, args_list)
+
+    freq = [1, 5] * apu.GHz
+    dist = np.array([5, 20, 1000, 20000])[:, np.newaxis] * apu.m
+
+    PL_los, PL_nlos, los_prob = _func(freq, dist)
+    print(PL_los)
+    print(PL_nlos)
+    print(los_prob)
+    assert_quantity_allclose(
+        PL_los,
+        [
+            [      np.nan,       np.nan],
+            [ 60.76626079,  74.74566088],
+            [108.24721393, 109.7252045 ],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )
+    assert_quantity_allclose(
+        PL_nlos,
+        [
+            [      np.nan,       np.nan],
+            [ 71.74479418,  85.72419426],
+            [130.78468516, 144.76408525],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )
+    assert_quantity_allclose(
+        los_prob,
+        [
+            [1.00000000e+00, 1.00000000e+00],
+            [9.72799557e-01, 9.72799557e-01],
+            [1.80001255e-02, 1.80001255e-02],
+            [9.00000000e-04, 9.00000000e-04],
+            ] * cnv.dimless,
+        )
+
+    PL_los, PL_nlos, los_prob = _func(freq, dist, h_bs=20 * apu.m)
+    print(PL_los)
+    print(PL_nlos)
+    assert_quantity_allclose(
+        PL_los,
+        [
+            [      np.nan,       np.nan],
+            [ 59.57605227,  73.55545236],
+            [110.07255004, 111.54965644],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )
+    assert_quantity_allclose(
+        PL_nlos,
+        [
+            [      np.nan,       np.nan],
+            [ 69.63055103,  83.60995112],
+            [130.78290388, 144.76230396],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )
+
+
+def test_imt_urban_micro_losses():
+
+    _func = pathprof.imt_urban_micro_losses
+    args_list = [
+        (0.5, 100, apu.GHz),
+        (0, 100000, apu.m),
+        (10, 150, apu.m),
+        (1, 22.5, apu.m),
+        ]
+
+    check_astro_quantities(_func, args_list)
+
+    freq = [1, 5] * apu.GHz
+    dist = np.array([5, 20, 1000, 20000])[:, np.newaxis] * apu.m
+
+    PL_los, PL_nlos, los_prob = _func(freq, dist)
+    print(PL_los)
+    print(PL_nlos)
+    assert_quantity_allclose(
+        PL_los,
+        [
+            [      np.nan,       np.nan],
+            [ 60.47880565,  74.45820574],
+            [118.53377126, 119.31141301],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )
+    assert_quantity_allclose(
+        PL_nlos,
+        [
+            [      np.nan,       np.nan],
+            [ 69.59913521,  84.4871963 ],
+            [128.3005538 , 143.18861489],
+            [      np.nan,       np.nan],
+            ] * cnv.dB,
+        )