diff --git a/pycraf/pathprof/cyprop.pyx b/pycraf/pathprof/cyprop.pyx
index 85a299e8..c09e14d9 100644
--- a/pycraf/pathprof/cyprop.pyx
+++ b/pycraf/pathprof/cyprop.pyx
@@ -150,7 +150,7 @@ PARAMETERS_BASIC = [
     ]
 
 
-PARAMETERS_V16 = [
+PARAMETERS_V18 = PARAMETERS_V16 = [
     ('path_type_50', '12d', '(0 - LOS, 1 - transhoriz)', cnv.dimless),
     ('d_bp_50', '12.6f', 'km', apu.km),
     ('h_bp_50', '12.6f', 'm', apu.m),
@@ -191,6 +191,9 @@ PARAMETERS_V16 = [
     ('S_tr_zh_b0', '12.6f', 'm / km', apu.m / apu.km),
     ]
 
+PARAMETERS_V18.append(('hS_tim_50', '12.6f', 'm / km', apu.m / apu.km))
+PARAMETERS_V18.append(('hS_tim_b0', '12.6f', 'm / km', apu.m / apu.km))
+
 
 PARAMETERS_V14 = [
     ('zeta_m', '12.6f', 'dimless', cnv.dimless),
@@ -286,6 +289,7 @@ cdef struct ppstruct:
     double nu_bull_50  # dimless
     int nu_bull_idx_50  # dimless
     double S_tim_50  # m / km
+    double hS_tim_50  # m / km  # version 18 only
     double S_rim_50  # m / km
     double S_tr_50  # m / km
 
@@ -296,6 +300,7 @@ cdef struct ppstruct:
     double nu_bull_b0  # dimless
     int nu_bull_idx_b0  # dimless
     double S_tim_b0  # m / km
+    double hS_tim_b0  # m / km  # version 18 only
     double S_rim_b0  # m / km
     double S_tr_b0  # m / km
     # double a_e_zh_50  # km
@@ -398,8 +403,11 @@ cdef class _PathProp(object):
             # override if you don't want builtin method:
             delta_N, N0,
             # override if you don't want builtin method:
-            hprof_dists, hprof_heights,
-            hprof_bearing, hprof_backbearing,
+            # hprof_dists, hprof_heights,
+            # hprof_bearing, hprof_backbearing,
+            hprof_data=None,
+            # allow to define a clutter height function for version 18
+            clutter_height_func=None,
             # set terrain heights to zero if desired
             # (only if hprof_xxx set to None aka automatic)
             bint generic_heights=False,
@@ -407,7 +415,12 @@ cdef class _PathProp(object):
             ):
 
         assert time_percent <= 50.
-        assert version == 14 or version == 16
+        assert version == 14 or version == 16 or version == 18
+        if version == 18 and hprof_data is None:
+            assert clutter_height_func is not None, (
+                'for version 18, a "clutter_height_func" must be provided, '
+                'alternatively, "hprof_data" can be used with "cl_heights"'
+            )
 
         assert zone_t >= -1 and zone_t <= 11
         assert zone_r >= -1 and zone_r <= 11
@@ -416,13 +429,18 @@ cdef class _PathProp(object):
             'delta_N and N0 must both be None or both be provided'
             )
 
-        assert (
-            (hprof_dists is None) == (hprof_heights is None) ==
-            (hprof_bearing is None) == (hprof_backbearing is None)
-            ), (
-                'hprof_dists, hprof_heights, bearing, and back_bearing '
-                'must all be None or all be provided'
-                )
+        if hprof_data is not None:
+            assert 'dists' in hprof_data, (
+                '"hprof_data" needs to include a "dists" array')
+            assert 'heights' in hprof_data, (
+                '"hprof_data" needs to include a "heights" array')
+            assert 'bearing' in hprof_data, (
+                '"hprof_data" needs to include a "bearing" array')
+            assert 'backbearing' in hprof_data, (
+                '"hprof_data" needs to include a "backbearing" array')
+            if version == 18:
+                assert 'cl_heights' in hprof_data, (
+                    '"hprof_data" needs to include a "cl_heights" array')
 
         self._pp.version = version
         self._pp.freq = freq
@@ -450,7 +468,7 @@ cdef class _PathProp(object):
         self._pp.time_percent = time_percent
         self._pp.polarization = polarization
 
-        if hprof_dists is None:
+        if hprof_data is None:
             (
                 lons,
                 lats,
@@ -466,13 +484,19 @@ cdef class _PathProp(object):
                     hprof_step,
                     generic_heights=generic_heights
                     )
+            if version == 18 and not generic_heights:
+                clheights = clutter_height_func(lons, lats)
+            else:
+                clheights = np.zeros_like(heights)
         else:
-            distances = hprof_dists.astype(np.float64, order='C', copy=False)
-            heights = hprof_heights.astype(np.float64, order='C', copy=False)
+            distances = hprof_data['dists'].astype(np.float64, order='C', copy=False)
+            heights = hprof_data['heights'].astype(np.float64, order='C', copy=False)
             hsize = distances.size
             distance = distances[hsize - 1]
-            bearing = hprof_bearing
-            back_bearing = hprof_backbearing
+            bearing = hprof_data['bearing']
+            back_bearing = hprof_data['backbearing']
+            if 'cl_heights' in hprof_data:
+                clheights = hprof_data['cl_heights']
 
         if len(distances) < 5:
             raise ValueError('Height profile must have at least 5 steps.')
@@ -509,7 +533,7 @@ cdef class _PathProp(object):
         hsize = distances.size
         mid_idx = hsize // 2
 
-        if hprof_dists is None:
+        if hprof_data is None:
             self._pp.lon_mid = lons[mid_idx]
             self._pp.lat_mid = lats[mid_idx]
         else:
@@ -540,6 +564,7 @@ cdef class _PathProp(object):
             distances,
             heights,
             zheights,
+            clheights,
             )
 
 
@@ -616,6 +641,7 @@ cdef void _process_path(
         double[::1] distances_view,
         double[::1] heights_view,
         double[::1] zheights_view,
+        double[::1] clheights_view,
         # double bearing,
         # double back_bearing,
         # double distance,
@@ -669,6 +695,29 @@ cdef void _process_path(
     pp.a_e_50 = 6371. * 157. / (157. - pp.delta_N)
     pp.a_e_b0 = 6371. * 3.
 
+    if pp.version == 18:
+        (
+            pp.path_type_50, pp.d_bp_50, pp.h_bp_50, pp.h_eff_50,
+            pp.nu_bull_50, pp.nu_bull_idx_50,
+            pp.S_tim_50, pp.hS_tim_50, pp.S_rim_50, pp.S_tr_50
+            ) = _diffraction_helper_v18(
+            pp.a_e_50, pp.distance,
+            distances_view, heights_view, clheights_view,
+            pp.h_ts, pp.h_rs,
+            pp.wavelen,
+            )
+
+        (
+            pp.path_type_b0, pp.d_bp_b0, pp.h_bp_b0, pp.h_eff_b0,
+            pp.nu_bull_b0, pp.nu_bull_idx_b0,
+            pp.S_tim_b0, pp.hS_tim_b0, pp.S_rim_b0, pp.S_tr_b0
+            ) = _diffraction_helper_v18(
+            pp.a_e_b0, pp.distance,
+            distances_view, heights_view, clheights_view,
+            pp.h_ts, pp.h_rs,
+            pp.wavelen,
+            )
+
     if pp.version == 16:
         (
             pp.path_type_50, pp.d_bp_50, pp.h_bp_50, pp.h_eff_50,
@@ -693,6 +742,7 @@ cdef void _process_path(
             )
 
         # similarly, we have to repeat the game with heights set to zero
+    if pp.version == 16 or pp.version == 18:
 
         (
             pp.path_type_zh_50, pp.d_bp_zh_50, pp.h_bp_zh_50, pp.h_eff_zh_50,
@@ -744,7 +794,7 @@ cdef void _process_path(
 
     if pp.version == 14:
         diff_edge_idx = pp.i_m50
-    elif pp.version == 16:
+    elif pp.version == 16 or pp.version == 18:
         diff_edge_idx = pp.nu_bull_idx_50
 
     (
@@ -853,6 +903,125 @@ cdef (double, double) _effective_antenna_heights(
     return (h_std, h_srd)
 
 
+cdef (
+    int, double, double, double, double, int, double, double, double, double
+    ) _diffraction_helper_v18(
+        double a_p,
+        double distance,
+        double[::1] d_v,
+        double[::1] h_v,
+        double[::1] c_v,  # clutter heights
+        double h_ts, double h_rs,
+        double wavelen,
+        ) noexcept nogil:
+
+    cdef:
+        int i, dsize
+        double d = distance, lam = wavelen, C_e500 = 500. / a_p
+        int path_type
+
+        # heights with clutter ("g[i]" in P.452-18)
+        double slope_i, slope_j, S_tim = -1.e31, S_tr, S_rim = -1.e31
+        # heights without clutter ("h[i]" in P.452-18)
+        double hslope_i, hS_tim = -1.e31
+
+        int nu_bull_idx
+        double d_bp, nu_bull = -1.e31, nu_i
+        double h_bp, h_eff, h_eff_i
+        double x, y  # temporary vars
+
+    dsize = d_v.shape[0]
+
+    for i in range(1, dsize - 1):
+
+        hslope_i = (
+            h_v[i] + C_e500 * d_v[i] * (d - d_v[i]) - h_ts
+            ) / d_v[i]
+        if (d_v[i] - d_v[0] > 0.05) and (d_v[dsize - 1] - d_v[i] > 0.05):
+            slope_i = (
+                h_v[i] + c_v[i] + C_e500 * d_v[i] * (d - d_v[i]) - h_ts
+                ) / d_v[i]
+        else:
+            slope_i = hslope_i
+
+
+        if slope_i > S_tim:
+            S_tim = slope_i
+
+        if hslope_i > hS_tim:
+            hS_tim = hslope_i
+
+    S_tr = (h_rs - h_ts) / d
+
+    if S_tim < S_tr:
+        path_type = 0
+    else:
+        path_type = 1
+
+    if path_type == 1:
+        # transhorizon
+        # find Bullington point, etc.
+        for i in range(1, dsize - 1):
+
+            if (d_v[i] - d_v[0] > 0.05) and (d_v[dsize - 1] - d_v[i] > 0.05):
+                slope_j = (
+                    h_v[i] + c_v[i] + C_e500 * d_v[i] * (d - d_v[i]) - h_rs
+                    ) / (d - d_v[i])
+            else:
+                slope_j = (
+                    h_v[i] + C_e500 * d_v[i] * (d - d_v[i]) - h_rs
+                    ) / (d - d_v[i])
+
+            if slope_j > S_rim:
+                S_rim = slope_j
+
+        d_bp = x = (h_rs - h_ts + S_rim * d) / (S_tim + S_rim)
+        y = a_p + h_ts / 1000 + d_bp * (S_tim / 1000 - d / 2 / a_p)
+        h_bp = 1000 * (sqrt(x ** 2 + y ** 2) - a_p)
+
+        h_eff = (
+            h_ts + S_tim * d_bp -
+            (
+                h_ts * (d - d_bp) + h_rs * d_bp
+                ) / d
+            )
+
+        nu_bull = h_eff * sqrt(
+            0.002 * d / lam / d_bp / (d - d_bp)
+            )  # == nu_b in Eq. 20
+        nu_bull_idx = -1  # dummy value
+
+    else:
+        # LOS
+
+        # find Bullington point, etc.
+
+        S_rim = NAN
+
+        # diffraction parameter
+        for i in range(1, dsize - 1):
+            h_eff_i = (
+                h_v[i] +
+                C_e500 * d_v[i] * (d - d_v[i]) -
+                (h_ts * (d - d_v[i]) + h_rs * d_v[i]) / d
+                )
+            nu_i = h_eff_i * sqrt(
+                0.002 * d / lam / d_v[i] / (d - d_v[i])
+                )
+            if nu_i > nu_bull:
+                nu_bull = nu_i
+                nu_bull_idx = i
+                h_eff = h_eff_i
+
+        d_bp = x = d_v[nu_bull_idx]
+        y = a_p + h_ts / 1000 + d_bp * (S_tr / 1000 - d / 2 / a_p)
+        h_bp = 1000 * (sqrt(x ** 2 + y ** 2) - a_p)
+
+    return (
+        path_type, d_bp, h_bp, h_eff, nu_bull, nu_bull_idx, S_tim, hS_tim, S_rim, S_tr
+        )
+
+
 cdef (
     int, double, double, double, double, int, double, double, double
     ) _diffraction_helper_v16(
@@ -1762,7 +1931,7 @@ cdef (double, double, double, double, double) _diffraction_loss_complete(
         double L_bfsg, E_sp, E_sbeta, L_min_b0p
         double F_i
 
-    if pp.version == 16:
+    if pp.version == 16 or pp.version == 18:
         L_d_50 = _delta_bullington_loss(pp, 0)
     elif pp.version == 14:
         L_d_50 = _diffraction_deygout_helper(
@@ -1780,7 +1949,7 @@ cdef (double, double, double, double, double) _diffraction_loss_complete(
 
     else:
 
-        if pp.version == 16:
+        if pp.version == 16 or pp.version == 18:
             L_d_beta = _delta_bullington_loss(pp, 1)
         elif pp.version == 14:
             L_d_beta = _diffraction_deygout_helper(
@@ -1868,6 +2037,7 @@ cdef (double, double, double, double, double, double, double) _path_attenuation_
         double L_min_bap, L_bam, L_b, L_b_corr
 
         double A_ht = 0., A_hr = 0.
+        double S_tim_tmp, S_tr_tmp
 
     if pp.zone_t != CLUTTER.UNKNOWN:
         A_ht = _clutter_correction(
@@ -1879,14 +2049,37 @@ cdef (double, double, double, double, double, double, double) _path_attenuation_
             )
 
     # not sure, if the 50% S_tim and S_tr values are to be used here...
-    if pp.version == 16:
+    if pp.version == 18:
+
+        if pp.time_percent > pp.beta0:
+            S_tim_tmp = pp.hS_tim_50
+            S_tr_tmp = pp.S_tr_50
+        else:
+            S_tim_tmp = pp.hS_tim_b0
+            S_tr_tmp = pp.S_tr_b0
+
+        F_j = 1 - 0.5 * (1. + tanh(
+            3. * _XI * (S_tim_tmp - S_tr_tmp) / _THETA
+            ))
+
+    elif pp.version == 16:
+
+        if pp.time_percent > pp.beta0:
+            S_tim_tmp = pp.S_tim_50
+            S_tr_tmp = pp.S_tr_50
+        else:
+            S_tim_tmp = pp.S_tim_b0
+            S_tr_tmp = pp.S_tr_b0
+
         F_j = 1 - 0.5 * (1. + tanh(
-            3. * _XI * (pp.S_tim_50 - pp.S_tr_50) / _THETA
+            3. * _XI * (S_tim_tmp - S_tr_tmp) / _THETA
             ))
+
     elif pp.version == 14:
         F_j = 1 - 0.5 * (1. + tanh(
             3. * _XI * (pp.theta - _THETA) / _THETA
             ))
+
     F_k = 1 - 0.5 * (1. + tanh(
         3. * _KAPPA * (pp.distance - _D_SW) / _D_SW
         ))
@@ -2101,6 +2294,7 @@ def height_map_data_cython(
         d_ct=None, d_cr=None,
         omega=None,
         int do_lonlat_profs=0,
+        clutter_height_func=None,
         ):
 
     '''
@@ -2141,6 +2335,15 @@ def height_map_data_cython(
     do_lonlat_profs : int, optional
         If True, also add `lons_profs` and `lats_profs` to output dict.
         (See below for further information.)
+    clutter_height_func : function, optional
+        If provided, the function is called with location data (lon and lat
+        of each pixel) in order to store effective clutter heights into
+        `hprof_data` for further use (it's also possible to add that
+        manually later). Must have the following signature
+
+            clutter_height_func(double lons_deg, double lats_deg) --> double
+
+        and support numpy broadcasting for lons and lats.
 
     Returns
     -------
@@ -2245,6 +2448,16 @@ def height_map_data_cython(
           Zero-valued array of the same length as `height_profs` for
           convenience.
 
+        - "cl_height_profs" : `~numpy.ndarray` 2D (float, (me, mh))
+
+          Clutter heights array of the same dimension as `height_profs`.
+          Only stored if `clutter_height_func` was provided as an argument.
+          This is required to perform P.452-18 calculations (where clutter
+          heights are added to the terrain heights).
+
+          This array can also be added manually before calling
+          `atten_map_fast`.
+
         - "lons_profs" : `~numpy.ndarray` 2D (float, (me, mh))
 
           Longitude (profiles) to each of the pixels on the map edge,
@@ -2289,7 +2502,7 @@ def height_map_data_cython(
         np.float64_t[:, ::1] _dist_map
         np.float64_t[:, ::1] _bearing_map, _backbearing_map
 
-    print('using hprof_step = {:.1f} m'.format(hprof_step))
+    # print('using hprof_step = {:.1f} m'.format(hprof_step))
 
     cosdelta = 1. / cos(DEG2RAD * lat_t) if do_cos_delta else 1.
 
@@ -2307,10 +2520,10 @@ def height_map_data_cython(
     lon_t_rad, lat_t_rad = DEG2RAD * lon_t, DEG2RAD * lat_t
     xcoords_rad = np.radians(xcoords)
     ycoords_rad = np.radians(ycoords)
-    print(
-        xcoords[0], xcoords[len(xcoords) - 1],
-        ycoords[0], ycoords[len(ycoords) - 1]
-        )
+    # print(
+    #     xcoords[0], xcoords[len(xcoords) - 1],
+    #     ycoords[0], ycoords[len(ycoords) - 1]
+    #     )
 
     # find max distance (will be one of the edges)
     max_distance = max([
@@ -2343,7 +2556,7 @@ def height_map_data_cython(
             (xcoords_rad.size - 1, xcoords_rad.size - 2, ycoords_rad.size - 1)
             ]
         ]) / 2  # rad
-    print('min pos angle resolution (at corner coords)', min_pa_res)
+    # print('min pos angle resolution (at corner coords)', min_pa_res)
 
 
     # path_idx_map stores the index of the edge-path that is closest
@@ -2394,7 +2607,7 @@ def height_map_data_cython(
         0, max_distance + hprof_step, hprof_step
         )
 
-    print('get geoid positions and bearings')
+    # print('get geoid positions and bearings')
     lons_rad, lats_rad, back_bearings_rad = cygeodesics.direct_cython(
         lon_t_rad, lat_t_rad,
         start_bearings[:, np.newaxis],
@@ -2406,7 +2619,7 @@ def height_map_data_cython(
     # print(xcoords, ycoords)
     # print(lons.min(), lons.max(), lats.min(), lats.max())
 
-    print('srtm query')
+    # print('srtm query')
     # hgt_res may not be set correctly yet, if no call to srtm was made before
     # let's do a simple query to make sure, it is set
     srtm._srtm_height_data(lon_t, lat_t)
@@ -2442,7 +2655,7 @@ def height_map_data_cython(
 
     refx, refy = _xcoords[0], _ycoords[0]
 
-    print('nogil loops')
+    # print('nogil loops')
 
     with nogil:
         for bidx in range(_start_bearings.shape[0]):
@@ -2558,6 +2771,9 @@ def height_map_data_cython(
     hprof_data['height_profs'] = height_profs
     hprof_data['zheight_prof'] = zheight_prof
 
+    if clutter_height_func is not None:
+        hprof_data['cl_height_profs'] = clutter_height_func(lons, lats)
+
     if do_lonlat_profs:
         hprof_data['lons_profs'] = lons
         hprof_data['lats_profs'] = lats
@@ -2598,7 +2814,9 @@ def atten_map_fast_cython(
         Polarization (default: 0)
         Allowed values are: 0 - horizontal, 1 - vertical
     version : int, optional
-        ITU-R Rec. P.452 version. Allowed values are: 14, 16
+        ITU-R Rec. P.452 version. Allowed values are: 14, 16, 18
+        If version 18 is used, `hprof_data` must have clutter heights
+        (`cl_height_profs`)
     base_water_density : double, optional
         Set base water level density (default: 7.5 g / m^3)
         See Rec. ITU-R P.452, Eq. (9a)
@@ -2647,7 +2865,7 @@ def atten_map_fast_cython(
     # clutter zone type is possible for each of Tx and Rx
 
     assert time_percent <= 50.
-    assert version == 14 or version == 16
+    assert version == 14 or version == 16 or version == 18
 
     cdef:
         # must set gains to zero, because gain is direction dependent
@@ -2698,11 +2916,15 @@ def atten_map_fast_cython(
 
         double[::1] dist_prof_v = _cf(hprof_data['dist_prof'])
         double[:, ::1] height_profs_v = _cf(hprof_data['height_profs'])
+        double[:, ::1] cl_height_profs_v
         double[::1] zheight_prof_v = _cf(hprof_data['zheight_prof'])
 
     xlen = len(xcoords)
     ylen = len(ycoords)
 
+    if version == 18:
+        cl_height_profs_v = _cf(hprof_data['cl_height_profs'])
+
     with nogil, parallel():
 
         pp = <ppstruct *> malloc(sizeof(ppstruct))
@@ -2770,15 +2992,23 @@ def atten_map_fast_cython(
                 # heights_v = height_profs_v[eidx, 0:didx + 1]
                 # zheights_v = zheight_prof_v[0:didx + 1]
 
-                _process_path(
-                    pp,
-                    # dists_v,
-                    # heights_v,
-                    # zheights_v,
-                    dist_prof_v[0:didx + 1],
-                    height_profs_v[eidx, 0:didx + 1],
-                    zheight_prof_v[0:didx + 1],
-                    )
+                if pp.version == 18:
+                    _process_path(
+                        pp,
+                        dist_prof_v[0:didx + 1],
+                        height_profs_v[eidx, 0:didx + 1],
+                        zheight_prof_v[0:didx + 1],
+                        cl_height_profs_v[eidx, 0:didx + 1],
+                        )
+                else:
+                    _process_path(
+                        pp,
+                        dist_prof_v[0:didx + 1],
+                        height_profs_v[eidx, 0:didx + 1],
+                        zheight_prof_v[0:didx + 1],
+                        zheight_prof_v[0:didx + 1],  # dummy
+                        )
+
 
                 (
                     L_b0p, L_bd, L_bs, L_ba, L_b, L_b_corr, L_dummy
@@ -2836,7 +3066,9 @@ def atten_path_fast_cython(
         Polarization (default: 0)
         Allowed values are: 0 - horizontal, 1 - vertical
     version : int, optional
-        ITU-R Rec. P.452 version. Allowed values are: 14, 16
+        ITU-R Rec. P.452 version. Allowed values are: 14, 16, 18
+        If version 18 is used, `hprof_data` must have clutter heights
+        (`cl_heights`)
     base_water_density : double, optional
         Set base water level density (default: 7.5 g / m^3)
         See Rec. ITU-R P.452, Eq. (9a)
@@ -2884,7 +3116,7 @@ def atten_path_fast_cython(
     # clutter zone type is possible for each of Tx and Rx
 
     assert time_percent <= 50.
-    assert version == 14 or version == 16
+    assert version == 14 or version == 16 or version == 18
 
     cdef:
         # must set gains to zero, because gain is direction dependent
@@ -2899,6 +3131,7 @@ def atten_path_fast_cython(
 
         double[::1] distances_v = _cf(hprof_data['distances'])
         double[::1] heights_v = _cf(hprof_data['heights'])
+        double[::1] cl_heights_v
         double[::1] zheights_v = np.zeros_like(_cf(hprof_data['heights']))
         double[::1] omega_v = _cf(hprof_data['omega'])
         double[::1] d_tm_v = _cf(hprof_data['d_tm'])
@@ -2916,6 +3149,9 @@ def atten_path_fast_cython(
     float_res = np.zeros((10, max_path_length), dtype=np.float64)
     int_res = np.zeros((1, max_path_length), dtype=np.int32)
 
+    if version == 18:
+        cl_heights_v = _cf(hprof_data['cl_heights'])
+
     cdef:
         double[:, :] float_res_v = float_res
         int[:, :] int_res_v = int_res
@@ -2980,12 +3216,22 @@ def atten_path_fast_cython(
 
             pp.distance = distances_v[i]
 
-            _process_path(
-                pp,
-                distances_v[0:i + 1],
-                heights_v[0:i + 1],
-                zheights_v[0:i + 1],
-                )
+            if pp.version == 18:
+                _process_path(
+                    pp,
+                    distances_v[0:i + 1],
+                    heights_v[0:i + 1],
+                    zheights_v[0:i + 1],
+                    cl_heights_v[0:i + 1],
+                    )
+            else:
+                _process_path(
+                    pp,
+                    distances_v[0:i + 1],
+                    heights_v[0:i + 1],
+                    zheights_v[0:i + 1],
+                    zheights_v[0:i + 1],  # dummy
+                    )
 
             (
                 L_b0p, L_bd, L_bs, L_ba, L_b, L_b_corr, L_dummy
@@ -3028,9 +3274,11 @@ def losses_complete_cython(
         # override if you don't want builtin method:
         delta_N=None, N0=None,
         # override if you don't want builtin method:
-        hprof_dists=None,
-        hprof_heights=None,
-        hprof_bearing=None, hprof_backbearing=None,
+        # hprof_dists=None,
+        # hprof_heights=None,
+        # hprof_bearing=None, hprof_backbearing=None,
+        hprof_data=None,
+        clutter_height_func=None,
         generic_heights=False,
         double base_water_density=7.5,  # g/m^3
         ):
@@ -3050,8 +3298,8 @@ def losses_complete_cython(
 
         # other
 
-        np.ndarray[double] lons, lats, distances, heights, zheights
-        double[::1] distances_v, heights_v, zheights_v
+        np.ndarray[double] lons, lats, distances, heights, zheights, clheights
+        double[::1] distances_v, heights_v, zheights_v, clheights_v
         double distance, bearing, back_bearing
         double lon_mid, lat_mid
         double _delta_N, _N0
@@ -3074,7 +3322,12 @@ def losses_complete_cython(
         int i, size
 
     assert np.all(time_percent <= 50.)
-    assert np.all((version == 14) | (version == 16))
+    assert np.all((version == 14) | (version == 16) | (version == 18))
+    if np.any(version == 18) and hprof_data is None:
+        assert clutter_height_func is not None, (
+            'for version 18, a "clutter_height_func" must be provided, '
+            'alternatively, "hprof_data" can be used with "cl_heights"'
+        )
 
     assert np.all((zone_t >= -1) & (zone_t <= 11))
     assert np.all((zone_r >= -1) & (zone_r <= 11))
@@ -3083,15 +3336,20 @@ def losses_complete_cython(
         'delta_N and N0 must both be None or both be provided'
         )
 
-    assert (
-        (hprof_dists is None) == (hprof_heights is None) ==
-        (hprof_bearing is None) == (hprof_backbearing is None)
-        ), (
-            'hprof_dists, hprof_heights, bearing, and back_bearing '
-            'must all be None or all be provided'
-            )
-
-    if hprof_dists is None:
+    if hprof_data is not None:
+        assert 'dists' in hprof_data, (
+            '"hprof_data" needs to include a "dists" array')
+        assert 'heights' in hprof_data, (
+            '"hprof_data" needs to include a "heights" array')
+        assert 'bearing' in hprof_data, (
+            '"hprof_data" needs to include a "bearing" array')
+        assert 'backbearing' in hprof_data, (
+            '"hprof_data" needs to include a "backbearing" array')
+        if np.any(version == 18):
+            assert 'cl_heights' in hprof_data, (
+                '"hprof_data" needs to include a "cl_heights" array')
+
+    if hprof_data is None:
         (
             lons, lats, distance, distances, heights,
             bearing, back_bearing, _,
@@ -3101,13 +3359,19 @@ def losses_complete_cython(
                 hprof_step,
                 generic_heights=generic_heights,
                 )
+        if np.any(version == 18) and not generic_heights:
+            clheights = clutter_height_func(lons, lats)
+        else:
+            clheights = np.zeros_like(heights)
     else:
-        distances = hprof_dists.astype(np.float64, order='C', copy=False)
-        heights = hprof_heights.astype(np.float64, order='C', copy=False)
+        distances = hprof_data['dists'].astype(np.float64, order='C', copy=False)
+        heights = hprof_data['heights'].astype(np.float64, order='C', copy=False)
         hsize = distances.size
         distance = distances[hsize - 1]
-        bearing = hprof_bearing
-        back_bearing = hprof_backbearing
+        bearing = hprof_data['bearing']
+        back_bearing = hprof_data['backbearing']
+        if 'cl_heights' in hprof_data:
+            clheights = hprof_data['cl_heights']
 
     if len(distances) < 5:
         raise ValueError('Height profile must have at least 5 steps.')
@@ -3116,12 +3380,13 @@ def losses_complete_cython(
 
     distances_v = distances
     heights_v = heights
+    clheights_v = clheights
     zheights_v = zheights
 
     hsize = distances.size
     mid_idx = hsize // 2
 
-    if hprof_dists is None:
+    if hprof_data is None:
         lon_mid = lons[mid_idx]
         lat_mid = lats[mid_idx]
     else:
@@ -3292,6 +3557,7 @@ def losses_complete_cython(
                         distances_v,
                         heights_v,
                         zheights_v,
+                        clheights_v,
                         )
 
                 pp.temperature = _temp[i]
diff --git a/pycraf/pathprof/propagation.py b/pycraf/pathprof/propagation.py
index 9919abf4..30870686 100644
--- a/pycraf/pathprof/propagation.py
+++ b/pycraf/pathprof/propagation.py
@@ -36,7 +36,7 @@
 # This is a wrapper class to expose the ppstruct members as attributes
 # (unfortunately, one cannot do dynamical attributes on cdef-classes)
 class PathProp(cyprop._PathProp):
-    '''
+    """
     Container class that holds all path profile properties.
 
     Parameters
@@ -78,7 +78,8 @@ class PathProp(cyprop._PathProp):
         Polarization (default: 0)
         Allowed values are: 0 - horizontal, 1 - vertical
     version : int, optional
-        ITU-R Rec. P.452 version. Allowed values are: 14, 16
+        ITU-R Rec. P.452 version. Allowed values are: 14, 16, 18
+        (With version 18, user has to provide clutter heights for the path)
     delta_N : `~astropy.units.Quantity`, optional
         Average radio-refractive index lapse-rate through the lowest 1 km of
         the atmosphere [N-units/km = 1/km]
@@ -86,22 +87,37 @@ class PathProp(cyprop._PathProp):
     N_0 : `~astropy.units.Quantity`, optional
         Sea-level surface refractivity [N-units = dimless]
         (default: query `~pycraf.pathprof.deltaN_N0_from_map`)
-    hprof_dists : `~astropy.units.Quantity`, optional
-        Distance vector associated with the height profile `hprof_heights`.
-        (default: query `~pycraf.pathprof.srtm_height_profile`)
-    hprof_heights : `~astropy.units.Quantity`, optional
-        Terrain heights profile for the distances in `hprof_dists`.
-        (default: query `~pycraf.pathprof.srtm_height_profile`)
-    hprof_bearing : `~astropy.units.Quantity`, optional
-        Start bearing of the height profile path.
-        (default: query `~pycraf.pathprof.srtm_height_profile`)
-    hprof_backbearing : `~astropy.units.Quantity`, optional
-        Back-bearing of the height profile path.
-        (default: query `~pycraf.pathprof.srtm_height_profile`)
+    hprof_data : dict, optional
+        This dictionary can be used to provide path profile data manually,
+        which can be more convenient (and in some cases faster) than
+        letting `PathProp` querying those data during instantiation
+        (with `~pycraf.pathprof.srtm_height_profile`). If provided, the
+        `hprof_data` takes precedence over other means to query profile data.
+        The dictionary must have the following entries:
+        - `dists` : `~numpy.array`
+            Distance vector associated with the height profile. [km]
+        - `heights` : `~numpy.array`
+            Terrain heights profile for the distances in `dists`. [m]
+        - `cl_heights` : `~numpy.array`, optional
+            Effective clutter heights associated with `heights` array. [m]
+            This is only required, if version 18 of P.452 is desired.
+        - `bearing` : double
+            Start bearing of the height profile path. [deg]
+        - `backbearing` : double
+            Back-bearing of the height profile path. [deg]
+    clutter_height_func : function, optional
+        If provided, the function is called with location data (lon and lat
+        of each pixel) in order to store effective clutter heights into
+        `hprof_data` for further use (it's also possible to add that
+        manually later). Must have the following signature
+
+            clutter_height_func(double lons_deg, double lats_deg) --> double
+
+        and support numpy broadcasting for lons and lats.
     generic_heights : bool
         If `generic_heights` is set to True, heights will be set to zero.
         This can be useful for generic (aka flat-Earth) computations.
-        The option is only meaningful, if the hprof_xxx parameters are set
+        The option is only meaningful, if `hprof_data` is set
         to `None` (which means automatic querying of the profiles).
         (Default: False)
     base_water_density : `~astropy.units.Quantity`, optional
@@ -117,7 +133,10 @@ class PathProp(cyprop._PathProp):
     -----
     - The diffraction-loss algorithm was changed between ITU-R P.452
       version 14 and 15. The former used a Deygout method, the new one
-      is based on a Bullington calculation with correction terms.
+      is based on a Bullington calculation with correction terms. Since
+      version 18, clutter heights are added to the terrain heights. The
+      user has to provide such data in the form of a function that
+      returns the requested heights for longitudes and latitudes.
 
     - Set `d_ct` and `d_cr` to zero for a terminal on ship or on a sea
       platform; only relevant if less than 5 km.
@@ -153,7 +172,7 @@ class PathProp(cyprop._PathProp):
       additional features such as automatic downloading of missing
       tiles or applying different interpolation methods (e.g., splines).
       For details see :ref:`working_with_srtm`.
-    '''
+    """
 
     @utils.ranged_quantity_input(
         freq=(0.1, 100, apu.GHz),
@@ -174,10 +193,10 @@ class PathProp(cyprop._PathProp):
         d_cr=(None, None, apu.m),
         delta_N=(None, None, cnv.dimless / apu.km),
         N0=(None, None, cnv.dimless),
-        hprof_dists=(None, None, apu.km),
-        hprof_heights=(None, None, apu.m),
-        hprof_bearing=(None, None, apu.deg),
-        hprof_backbearing=(None, None, apu.deg),
+        # hprof_dists=(None, None, apu.km),
+        # hprof_heights=(None, None, apu.m),
+        # hprof_bearing=(None, None, apu.deg),
+        # hprof_backbearing=(None, None, apu.deg),
         base_water_density=(0.5, 100, apu.g / apu.m ** 3),
         strip_input_units=True, allow_none=True, output_unit=None
         )
@@ -200,8 +219,10 @@ def __init__(
             # override if you don't want builtin method:
             delta_N=None, N0=None,
             # override if you don't want builtin method:
-            hprof_dists=None, hprof_heights=None,
-            hprof_bearing=None, hprof_backbearing=None,
+            # hprof_dists=None, hprof_heights=None,
+            # hprof_bearing=None, hprof_backbearing=None,
+            hprof_data=None,
+            clutter_height_func=None,
             generic_heights=False,
             base_water_density=7.5 * apu.g / apu.m ** 3
             ):
@@ -222,10 +243,12 @@ def __init__(
             polarization=polarization,
             version=version,
             delta_N=delta_N, N0=N0,
-            hprof_dists=hprof_dists,
-            hprof_heights=hprof_heights,
-            hprof_bearing=hprof_bearing,
-            hprof_backbearing=hprof_backbearing,
+            # hprof_dists=hprof_dists,
+            # hprof_heights=hprof_heights,
+            # hprof_bearing=hprof_bearing,
+            # hprof_backbearing=hprof_backbearing,
+            hprof_data=hprof_data,
+            clutter_height_func=clutter_height_func,
             generic_heights=generic_heights,
             base_water_density=base_water_density,
             )
@@ -235,6 +258,8 @@ def __init__(
             self.__params += cyprop.PARAMETERS_V14
         elif self._pp['version'] == 16:
             self.__params += cyprop.PARAMETERS_V16
+        elif self._pp['version'] == 18:
+            self.__params += cyprop.PARAMETERS_V18
 
         # no need to set property, as readonly and immutable
         # can just copy to __dict__
@@ -503,10 +528,11 @@ def height_map_data(
         d_ct=None, d_cr=None,
         omega_percent=0 * apu.percent,
         do_lonlat_profs=False,
+        clutter_height_func=None,
         cache_path=None, clobber=False,
         ):
 
-    '''
+    """
     Calculate height profiles and auxillary maps needed for
     `~pycraf.pathprof.atten_map_fast`.
 
@@ -565,6 +591,15 @@ def height_map_data(
     do_lonlat_profs : bool, optional
         If True, also add `lons_profs` and `lats_profs` to output dict.
         (See below for further information. Default: False)
+    clutter_height_func : function, optional
+        If provided, the function is called with location data (lon and lat
+        of each pixel) in order to store effective clutter heights into
+        `hprof_data` for further use (it's also possible to add that
+        manually later). Must have the following signature
+
+            clutter_height_func(double lons_deg, double lats_deg) --> double
+
+        and support numpy broadcasting for lons and lats.
     cache_path : str, optional
         If set, the `joblib package
         <https://joblib.readthedocs.io/en/latest/>`_ is used to cache
@@ -675,6 +710,16 @@ def height_map_data(
           Zero-valued array of the same length as `height_profs` for
           convenience.
 
+        - "cl_height_profs" : `~numpy.ndarray` 2D (float, (me, mh))
+
+          Clutter heights array of the same dimension as `height_profs`.
+          Only stored if `clutter_height_func` was provided as an argument.
+          This is required to perform P.452-18 calculations (where clutter
+          heights are added to the terrain heights).
+
+          This array can also be added manually before calling
+          `atten_map_fast`.
+
         - "lons_profs" : `~numpy.ndarray` 2D (float, (me, mh))
 
           Longitude (profiles) to each of the pixels on the map edge,
@@ -699,7 +744,7 @@ def height_map_data(
       additional features such as automatic downloading of missing
       tiles or applying different interpolation methods (e.g., splines).
       For details see :ref:`working_with_srtm`.
-    '''
+    """
 
     args = lon_t, lat_t, map_size_lon, map_size_lat
     kwargs = dict(
@@ -710,6 +755,7 @@ def height_map_data(
         d_ct=d_ct, d_cr=d_cr,
         omega=omega_percent,
         do_lonlat_profs=do_lonlat_profs,
+        clutter_height_func=clutter_height_func,
         )
 
     joblib_available = True
@@ -774,7 +820,7 @@ def atten_map_fast(
         version=16,
         base_water_density=7.5 * apu.g / apu.m ** 3,
         ):
-    '''
+    """
     Calculate attenuation maps using a fast method.
 
     Parameters
@@ -797,7 +843,9 @@ def atten_map_fast(
         Polarization (default: 0)
         Allowed values are: 0 - horizontal, 1 - vertical
     version : int, optional
-        ITU-R Rec. P.452 version. Allowed values are: 14, 16
+        ITU-R Rec. P.452 version. Allowed values are: 14, 16, 18
+        If version 18 is used, `hprof_data` must have clutter heights
+        (`cl_height_profs`)
     base_water_density : `~astropy.units.Quantity`, optional
         For atmospheric attenuation, the water vapor content plays a role.
         In Rec. ITU-R P.452, Eq. (9a), the water content is variable (depending on the fraction of the path over the water). However, the base level is set to :math:`7.5~\\mathrm{g}/\\mathrm{m}^3`. For extraordinarily dry places, which are often used for radio astronomy, this value can be too high.
@@ -842,7 +890,7 @@ def atten_map_fast(
       is based on a Bullington calculation with correction terms.
     - In future versions, more entries may be added to the results
       dictionary.
-    '''
+    """
 
     float_res, int_res = cyprop.atten_map_fast_cython(
         freq,
@@ -884,9 +932,10 @@ def height_path_data(
         lon_r, lat_r,
         step,
         zone_t=cyprop.CLUTTER.UNKNOWN, zone_r=cyprop.CLUTTER.UNKNOWN,
+        clutter_height_func=None,
         ):
 
-    '''
+    """
     Calculate height profile auxillary data needed for
     `~pycraf.pathprof.atten_path_fast`.
 
@@ -908,6 +957,15 @@ def height_path_data(
     zone_t, zone_r : CLUTTER enum, optional
         Clutter type for transmitter/receiver terminal.
         (default: CLUTTER.UNKNOWN)
+    clutter_height_func : function, optional
+        If provided, the function is called with location data (lon and lat
+        of each pixel) in order to store effective clutter heights into
+        `hprof_data` for further use (it's also possible to add that
+        manually later). Must have the following signature
+
+            clutter_height_func(double lons_deg, double lats_deg) --> double
+
+        and support numpy broadcasting for lons and lats.
 
     Returns
     -------
@@ -941,6 +999,16 @@ def height_path_data(
 
           Height profile.
 
+
+        - "cl_heights" : `~numpy.ndarray` 2D (float, (me, mh))
+
+          Clutter heights array of the same dimension as `heights`.
+          Only stored if `clutter_height_func` was provided as an argument.
+          This is required to perform P.452-18 calculations (where clutter
+          heights are added to the terrain heights).
+
+          This array can also be added manually before calling
+          `atten_path_fast`.
         - "bearing" : float
 
           Start bearing of path.
@@ -986,7 +1054,7 @@ def height_path_data(
     different clutter types (in the array). You can modify the returned
     arrays in `hprof_data`, of course, before feeding into
     `~pycraf.pathprof.atten_path_fast`.
-    '''
+    """
 
     (
         lons, lats, distance, distances, heights,
@@ -1034,6 +1102,10 @@ def height_path_data(
     hprof_data['N0'] = N0
     hprof_data['beta0'] = beta0
 
+    if clutter_height_func is not None:
+        hprof_data['cl_heights'] = clutter_height_func(lons, lats)
+
+
     return hprof_data
 
 
@@ -1189,6 +1261,9 @@ def height_path_data_generic(
     hprof_data['N0'] = N0
     hprof_data['beta0'] = beta0
 
+    # make sure that version 18 is working
+    hprof_data['cl_heights'] = np.zeros_like(heights)
+
     return hprof_data
 
 
@@ -1383,10 +1458,10 @@ def atten_path_fast(
     d_cr=(None, None, apu.m),
     delta_N=(None, None, cnv.dimless / apu.km),
     N0=(None, None, cnv.dimless),
-    hprof_dists=(None, None, apu.km),
-    hprof_heights=(None, None, apu.m),
-    hprof_bearing=(None, None, apu.deg),
-    hprof_backbearing=(None, None, apu.deg),
+    # hprof_dists=(None, None, apu.km),
+    # hprof_heights=(None, None, apu.m),
+    # hprof_bearing=(None, None, apu.deg),
+    # hprof_backbearing=(None, None, apu.deg),
     base_water_density=(0.5, 100, apu.g / apu.m ** 3),
     strip_input_units=True, allow_none=True, output_unit=None
     )
@@ -1409,12 +1484,14 @@ def losses_complete(
         # override if you don't want builtin method:
         delta_N=None, N0=None,
         # override if you don't want builtin method:
-        hprof_dists=None, hprof_heights=None,
-        hprof_bearing=None, hprof_backbearing=None,
+        # hprof_dists=None, hprof_heights=None,
+        # hprof_bearing=None, hprof_backbearing=None,
+        hprof_data=None,
+        clutter_height_func=None,
         generic_heights=False,
         base_water_density=7.5 * apu.g / apu.m ** 3,
         ):
-    '''
+    """
     Calculate propagation losses for a fixed path using a parallelized method.
 
     The difference to the usual `~pycraf.pathprof.PathProp` +
@@ -1475,18 +1552,33 @@ def losses_complete(
     N_0 : `~astropy.units.Quantity`, scalar, optional
         Sea-level surface refractivity [N-units = dimless]
         (default: query `~pycraf.pathprof.deltaN_N0_from_map`)
-    hprof_dists : `~astropy.units.Quantity`, optional
-        Distance vector associated with the height profile `hprof_heights`.
-        (default: query `~pycraf.pathprof.srtm_height_profile`)
-    hprof_heights : `~astropy.units.Quantity`, optional
-        Terrain heights profile for the distances in `hprof_dists`.
-        (default: query `~pycraf.pathprof.srtm_height_profile`)
-    hprof_bearing : `~astropy.units.Quantity`, optional
-        Start bearing of the height profile path.
-        (default: query `~pycraf.pathprof.srtm_height_profile`)
-    hprof_backbearing : `~astropy.units.Quantity`, optional
-        Back-bearing of the height profile path.
-        (default: query `~pycraf.pathprof.srtm_height_profile`)
+    hprof_data : dict, optional
+        This dictionary can be used to provide path profile data manually,
+        which can be more convenient (and in some cases faster) than
+        letting `PathProp` querying those data during instantiation
+        (with `~pycraf.pathprof.srtm_height_profile`). If provided, the
+        `hprof_data` takes precedence over other means to query profile data.
+        The dictionary must have the following entries:
+        - `dists` : `~numpy.array`
+            Distance vector associated with the height profile. [km]
+        - `heights` : `~numpy.array`
+            Terrain heights profile for the distances in `dists`. [m]
+        - `cl_heights` : `~numpy.array`, optional
+            Effective clutter heights associated with `heights` array. [m]
+            This is only required, if version 18 of P.452 is desired.
+        - `bearing` : double
+            Start bearing of the height profile path. [deg]
+        - `backbearing` : double
+            Back-bearing of the height profile path. [deg]
+    clutter_height_func : function, optional
+        If provided, the function is called with location data (lon and lat
+        of each pixel) in order to store effective clutter heights into
+        `hprof_data` for further use (it's also possible to add that
+        manually later). Must have the following signature
+
+            clutter_height_func(double lons_deg, double lats_deg) --> double
+
+        and support numpy broadcasting for lons and lats.
     generic_heights : bool
         If `generic_heights` is set to True, heights will be set to zero.
         This can be useful for generic (aka flat-Earth) computations.
@@ -1586,7 +1678,7 @@ def losses_complete(
     - In future versions, more entries may be added to the results
       dictionary.
 
-    '''
+    """
 
     res = cyprop.losses_complete_cython(
         freq,
@@ -1608,10 +1700,12 @@ def losses_complete(
         polarization=polarization,
         version=version,
         delta_N=delta_N, N0=N0,
-        hprof_dists=hprof_dists,
-        hprof_heights=hprof_heights,
-        hprof_bearing=hprof_bearing,
-        hprof_backbearing=hprof_backbearing,
+        # hprof_dists=hprof_dists,
+        # hprof_heights=hprof_heights,
+        # hprof_bearing=hprof_bearing,
+        # hprof_backbearing=hprof_backbearing,
+        hprof_data=hprof_data,
+        clutter_height_func=clutter_height_func,
         generic_heights=generic_heights,
         base_water_density=base_water_density,
         )
diff --git a/pycraf/pathprof/tests/data/cases.zip b/pycraf/pathprof/tests/data/cases.zip
index 2327ba66..09ca4b14 100644
Binary files a/pycraf/pathprof/tests/data/cases.zip and b/pycraf/pathprof/tests/data/cases.zip differ
diff --git a/pycraf/pathprof/tests/data/cases_generic.zip b/pycraf/pathprof/tests/data/cases_generic.zip
index 21530378..b2c7239c 100644
Binary files a/pycraf/pathprof/tests/data/cases_generic.zip and b/pycraf/pathprof/tests/data/cases_generic.zip differ
diff --git a/pycraf/pathprof/tests/data/fastmap.zip b/pycraf/pathprof/tests/data/fastmap.zip
index b02a5999..08af686b 100644
Binary files a/pycraf/pathprof/tests/data/fastmap.zip and b/pycraf/pathprof/tests/data/fastmap.zip differ
diff --git a/pycraf/pathprof/tests/test_propagation.py b/pycraf/pathprof/tests/test_propagation.py
index 65bc37d7..62fd938a 100644
--- a/pycraf/pathprof/tests/test_propagation.py
+++ b/pycraf/pathprof/tests/test_propagation.py
@@ -45,12 +45,25 @@
     ]
 
 
+
+def dummy_clutter_heights(lons, lats, do_zero=False):
+    with NumpyRNGContext(seed=0):
+        rnd_idxs = np.random.randint(4, size=lons.size)
+    return np.array([0., 10., 15., 20.])[rnd_idxs].reshape(lons.shape)
+
+
+def dummy_clutter_zheights(lons, lats):
+    return np.zeros(lons.shape)
+
 @pytest.mark.remote_data
 @pytest.mark.usefixtures('srtm_handler')
 class TestPropagation:
 
     def setup_method(self):
 
+        def _dummy_clutter_heights(lons, lats):
+            return np.zeros(lons.shape)
+
         # TODO: add further test cases
 
         print(pathprof.srtm.SrtmConf)
@@ -70,7 +83,7 @@ def setup_method(self):
         self.cases_h_tgs = [50., 200.]
         self.cases_h_rgs = [50., 200.]
         self.cases_time_percents = [2, 10, 50]
-        self.cases_versions = [14, 16]
+        self.cases_versions = [14, 16, 18]
         self.cases_G_ts = [0, 20]
         self.cases_G_rs = [0, 30]
 
@@ -92,7 +105,7 @@ def setup_method(self):
                 [0.1, 1., 10.],
                 [(50, 50), (200, 200)],
                 [2, 10, 50],
-                [14, 16],
+                [14, 16, 18],
                 [(0, 0)],
             ))
 
@@ -119,6 +132,7 @@ def setup_method(self):
                 self.hprof_step,
                 time_percent * apu.percent,
                 version=version,
+                clutter_height_func=dummy_clutter_heights,
                 )
 
             self.pprops.append(pprop)
@@ -204,6 +218,7 @@ def test_pathprop(self):
                 self.hprof_step,
                 time_percent * apu.percent,
                 version=version,
+                clutter_height_func=dummy_clutter_heights,
                 )
 
             with ZipFile(self.cases_zip_name) as myzip:
@@ -360,6 +375,7 @@ def test_complete_losses(self):
             G_r=G_rs * cnv.dBi,
             omega=self.omega,
             version=versions,
+            clutter_height_func=dummy_clutter_heights,
             )
 
         for tup in np.nditer([
@@ -392,6 +408,7 @@ def test_height_map_data_h5py(self, tmp_path_factory):
             6.5 * apu.deg, 50.5 * apu.deg,
             900 * apu.arcsec, 900 * apu.arcsec,
             map_resolution=30 * apu.arcsec,
+            clutter_height_func=dummy_clutter_heights,
             )
 
         # also testing reading/writing from hdf5
@@ -486,6 +503,7 @@ def test_fast_atten_map_h5py(self, tmp_path_factory):
             # "with ZipFile" below; results will be in the tmp-dir
             # and can be added manually to the zipfile
             # tfile = str(zipdir.join(fname))
+            # tfile = str(zipdir / fname)
             # print('writing temporary files to', zipdir)
             # with h5py.File(tfile, 'w') as h5f:
             #     for k, v in results.items():
@@ -518,6 +536,7 @@ def test_height_map_data_npz(self, tmp_path_factory):
             6.5 * apu.deg, 50.5 * apu.deg,
             900 * apu.arcsec, 900 * apu.arcsec,
             map_resolution=30 * apu.arcsec,
+            clutter_height_func=dummy_clutter_heights,
             )
 
         # also testing reading/writing from hdf5
@@ -609,7 +628,7 @@ def test_fast_atten_map_npz(self, tmp_path_factory):
             # Also, if you want to create all at once, comment-out the
             # "with ZipFile" below; results will be in the tmp-dir
             # and can be added manually to the zipfile
-            # tfile = str(zipdir.join(fname))
+            # tfile = str(zipdir / fname)
             # np.savez(tfile, **results)
 
             with ZipFile(self.fastmap_zip_name) as myzip:
@@ -667,6 +686,11 @@ def test_atten_path_fast(self):
             ) = pathprof.srtm_height_profile(
                 lon_t, lat_t, lon_r, lat_r, hprof_step
                 )
+        distances = distances.to_value(apu.km)
+        heights = heights.to_value(apu.m)
+        bearing = bearing.to_value(apu.deg)
+        back_bearings = back_bearings.to_value(apu.deg)
+        cl_heights = dummy_clutter_heights(lons, lats)
 
         atten_path = np.zeros((6, len(distances)), dtype=np.float64)
         eps_pt_path = np.zeros((len(distances)), dtype=np.float64)
@@ -674,8 +698,16 @@ def test_atten_path_fast(self):
         d_lt_path = np.zeros((len(distances)), dtype=np.float64)
         d_lr_path = np.zeros((len(distances)), dtype=np.float64)
 
+
         for idx in range(6, len(distances)):
 
+            hprof_data = {
+                "dists": distances[: idx + 1],
+                "heights": heights[: idx + 1],
+                "bearing": bearing,
+                "backbearing": back_bearings[idx],
+                "cl_heights": cl_heights[: idx + 1],
+                }
             pprop = pathprof.PathProp(
                 freq,
                 temperature, pressure,
@@ -685,10 +717,12 @@ def test_atten_path_fast(self):
                 hprof_step,
                 time_percent,
                 zone_t=zone_t, zone_r=zone_r,
-                hprof_dists=distances[:idx + 1],
-                hprof_heights=heights[:idx + 1],
-                hprof_bearing=bearing,
-                hprof_backbearing=back_bearings[idx],
+                # hprof_dists=distances[:idx + 1],
+                # hprof_heights=heights[:idx + 1],
+                # hprof_bearing=bearing,
+                # hprof_backbearing=back_bearings[idx],
+                hprof_data=hprof_data,
+                # clutter_height_func=dummy_clutter_heights,
                 # delta_N=hprof_data['delta_N'][idx] * cnv.dimless / apu.km,
                 # N0=hprof_data['N0'][idx] * cnv.dimless,
                 )
@@ -735,7 +769,7 @@ def setup_method(self):
         self.cases_h_tgs = [50., 200.]
         self.cases_h_rgs = [50., 200.]
         self.cases_time_percents = [2, 10, 50]
-        self.cases_versions = [14, 16]
+        self.cases_versions = [14, 16, 18]
         self.cases_G_ts = [0, 20]
         self.cases_G_rs = [0, 30]
 
@@ -757,7 +791,7 @@ def setup_method(self):
                 [0.1, 1., 10.],
                 [(50, 50), (200, 200)],
                 [2, 10, 50],
-                [14, 16],
+                [14, 16, 18],
                 [(0, 0)],
             ))
 
@@ -782,6 +816,7 @@ def setup_method(self):
                 self.hprof_step,
                 time_percent * apu.percent,
                 version=version,
+                clutter_height_func=dummy_clutter_zheights,
                 generic_heights=True,
                 )
 
@@ -868,6 +903,7 @@ def test_pathprop(self):
                 self.hprof_step,
                 time_percent * apu.percent,
                 version=version,
+                clutter_height_func=dummy_clutter_zheights,
                 generic_heights=True,
                 )
 
@@ -1069,6 +1105,11 @@ def test_atten_path_fast_generic(self):
                 lon_t, lat_mid, lon_r, lat_mid, hprof_step,
                 generic_heights=True,
                 )
+        distances = distances.to_value(apu.km)
+        heights = heights.to_value(apu.m)
+        bearing = bearing.to_value(apu.deg)
+        back_bearings = back_bearings.to_value(apu.deg)
+        cl_heights = np.zeros_like(heights)
 
         zone_t, zone_r = pathprof.CLUTTER.URBAN, pathprof.CLUTTER.SUBURBAN
 
@@ -1091,6 +1132,14 @@ def test_atten_path_fast_generic(self):
 
         for idx in range(6, len(distances)):
 
+            this_hprof_data = {
+                "dists": distances[: idx + 1],
+                "heights": heights[: idx + 1],
+                "bearing": bearing,
+                "backbearing": back_bearings[idx],
+                "cl_heights": cl_heights[: idx + 1],
+                }
+
             pprop = pathprof.PathProp(
                 freq,
                 temperature, pressure,
@@ -1100,10 +1149,12 @@ def test_atten_path_fast_generic(self):
                 hprof_step,
                 time_percent,
                 zone_t=zone_t, zone_r=zone_r,
-                hprof_dists=distances[:idx + 1],
-                hprof_heights=0 * heights[:idx + 1],
-                hprof_bearing=bearing,
-                hprof_backbearing=back_bearings[idx],
+                # hprof_dists=distances[:idx + 1],
+                # hprof_heights=0 * heights[:idx + 1],
+                # hprof_bearing=bearing,
+                # hprof_backbearing=back_bearings[idx],
+                hprof_data=this_hprof_data,
+                # clutter_height_func=dummy_clutter_zheights,
                 delta_N=hprof_data['delta_N'][idx] * cnv.dimless / apu.km,
                 N0=hprof_data['N0'][idx] * cnv.dimless,
                 )