From 97321ce898420633a4fe0d72a7117d84b35b071a Mon Sep 17 00:00:00 2001
From: iportillo <portillo@mit.edu>
Date: Wed, 30 May 2018 03:06:18 -0400
Subject: [PATCH] Merge for release 0.2.0 (#2)

* Update tests to include coverage of recommendation ITU-R P.618
* Reducing amount of code duplicated.
* Added recommendation ITU-R P.1853-0 and ITU-R P.1853-1
* Added Recommendation ITU-R P.530-17
* Added method to compute scintillation_attenuation_sigma
* Added method slant_inclined_path_equivalent_height
* Added documentation for plot_in_map function
* Added examples for recommendations ITU-1510 and ITU-837
---
 examples/itu_1510_mean_surface_temperature.py |  27 ++
 examples/itu_837_rainfall_map.py              |  28 ++
 itur/__init__.py                              |   2 +-
 itur/__version__.py                           |   2 +-
 itur/models/itu1853.py                        | 445 ++++++++++++++----
 itur/models/itu530.py                         | 137 ++++--
 itur/models/itu618.py                         |  94 +++-
 itur/models/itu676.py                         |  53 ++-
 itur/utils.py                                 |  44 +-
 setup.py                                      |   2 +-
 10 files changed, 698 insertions(+), 136 deletions(-)
 create mode 100644 examples/itu_1510_mean_surface_temperature.py
 create mode 100644 examples/itu_837_rainfall_map.py

diff --git a/examples/itu_1510_mean_surface_temperature.py b/examples/itu_1510_mean_surface_temperature.py
new file mode 100644
index 0000000..814cc0a
--- /dev/null
+++ b/examples/itu_1510_mean_surface_temperature.py
@@ -0,0 +1,27 @@
+# -*- coding: utf-8 -*-
+""" This example shows how to compute the anuaml mean surface temperature
+over a large region of the Earth.
+
+This image is similar to the one plotted in page 3 of Recommendation
+ITU-R P.1510-1.
+"""
+import itur
+import matplotlib.pyplot as plt
+
+# Set Recommendation ITU-R P.837 to version 7
+itur.models.itu1510.change_version(1)
+
+# Generate a regular grid of latitude and longitudes with 0.1 degree resolution
+# for the region of interest.
+lat, lon = itur.utils.regular_lat_lon_grid(resolution_lat=0.1,
+                                           resolution_lon=0.1)
+
+# Compute the surface mean temperature
+T = itur.models.itu1510.surface_mean_temperature(lat, lon)
+
+# Display the results in a map
+fig = plt.figure(figsize=(16, 8))
+ax = fig.add_subplot(1, 1, 1)
+m = itur.utils.plot_in_map(
+        T, lat, lon, cmap='jet', vmin=230, vmax=310, ax=ax,
+        cbar_text='Annual mean surface temperature [K]')
diff --git a/examples/itu_837_rainfall_map.py b/examples/itu_837_rainfall_map.py
new file mode 100644
index 0000000..b7958cd
--- /dev/null
+++ b/examples/itu_837_rainfall_map.py
@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+""" This example shows how to compute the rainfall-rate (mm/hr) exceeded
+for 0.01 % of the time of the average year over a large region of the Earth.
+
+This image is similar to the one plotted in page 5 of Recommendation
+ITU-R P.837-7.
+"""
+import itur
+import matplotlib.pyplot as plt
+
+# Set Recommendation ITU-R P.837 to version 7
+itur.models.itu837.change_version(7)
+
+# Generate a regular grid of latitude and longitudes with 0.1 degree resolution
+# for the region of interest.
+lat, lon = itur.utils.regular_lat_lon_grid(resolution_lat=0.1,
+                                           resolution_lon=0.1)
+
+# Compute the rainfall rate exceeded for 0.01 % of the time.
+p = 0.01
+R001 = itur.models.itu837.rainfall_rate(lat, lon, p)
+
+# Display the results in a map
+fig = plt.figure(figsize=(16, 8))
+ax = fig.add_subplot(1, 1, 1)
+m = itur.utils.plot_in_map(
+    R001, lat, lon, cmap='jet', vmin=0, vmax=90, ax=ax,
+    cbar_text='Rainfall rate exceeded for 0.01% of an average year [mm/hr]')
diff --git a/itur/__init__.py b/itur/__init__.py
index ca9da60..e9ab200 100644
--- a/itur/__init__.py
+++ b/itur/__init__.py
@@ -121,7 +121,7 @@ def atmospheric_attenuation_slant_path(
         Total atmospheric attenuation (dB)
 
     Ag, Ac, Ar, As, A : tuple
-        Gaseous, Cloud, Rain, Scintillation cotributions to total attenuation,
+        Gaseous, Cloud, Rain, Scintillation contributions to total attenuation,
         and total attenuation (dB)
 
 
diff --git a/itur/__version__.py b/itur/__version__.py
index 35017fa..1206a6f 100644
--- a/itur/__version__.py
+++ b/itur/__version__.py
@@ -2,4 +2,4 @@
 """Version information."""
 
 # The following line *must* be the last in the module, exactly as formatted:
-__version__ = "0.1.6"
+__version__ = "0.2.0"
diff --git a/itur/models/itu1853.py b/itur/models/itu1853.py
index 631671e..2f4f7e6 100644
--- a/itur/models/itu1853.py
+++ b/itur/models/itu1853.py
@@ -7,14 +7,18 @@
 import scipy.stats as stats
 from scipy.signal import lfilter
 
-from itu618 import fit_rain_attenuation_to_lognormal
-from itu840 import lognormal_approximation_coefficient, specific_attenuation_coefficients
+from itu618 import rain_attenuation, scintillation_attenuation_sigma
+from itu676 import gamma0_exact, \
+    slant_inclined_path_equivalent_height
+from itu840 import lognormal_approximation_coefficient, \
+    specific_attenuation_coefficients
+from itu835 import standard_pressure, standard_water_vapour_density
 from itu836 import total_water_vapour_content
-from itu837 import rain_percentage_probability
+from itu837 import rainfall_probability
 from itu676 import zenit_water_vapour_attenuation
 from itu1510 import surface_mean_temperature
-from itur.utils import load_data, dataset_dir, prepare_input_array,\
-    prepare_output_array, prepare_quantity
+from itu1511 import topographic_altitude
+from itur.utils import prepare_quantity
 
 from astropy import units as u
 
@@ -23,12 +27,8 @@ class __ITU1853():
     """Tropospheric attenuation time series synthesis
 
     Available versions include:
-    * P.1853-1 (08/94) (Superseded)
-    * P.1853-2 (10/99) (Superseded)
-    * P.1853-3 (02/01) (Superseded)
-    * P.1853-4 (04/03) (Superseded)
-    * P.1853-5 (08/07) (Superseded)
-    * P.1853-6 (02/12) (Current version)
+    * P.1853-0 (10/09) (Superseded)
+    * P.1853-1 (02/12) (Current version)
     """
     # This is an abstract class that contains an instance to a version of the
     # ITU-R P.1853 recommendation.
@@ -40,18 +40,40 @@ def __init__(self, version=1):
             self.instance = _ITU1853_0()
         else:
             raise ValueError(
-                'Version ' +
-                str(version) +
-                ' is not implemented' +
-                ' for the ITU-R P.1853 model.')
+                'Version {0} is not implemented for the ITU-R P.1853 model.'
+                .format(version))
+
+    def set_seed(self, seed):
+        np.random.seed(seed)
 
     @property
     def __version__(self):
         return self.instance.__version__
 
-    def rain_attenuation_synthesis(self, lat, lon, f, el, hs, T, Ts=1, n=None):
-        return self.instance.rain_attenuation_synthesis(lat, lon, f,
-                                                        el, hs, T, Ts=Ts, n=n)
+    def rain_attenuation_synthesis(self, lat, lon, f, el, hs, Ns,
+                                   Ts=1, n=None):
+        return self.instance.rain_attenuation_synthesis(
+                lat, lon, f, el, hs, Ns, Ts=Ts, n=n)
+
+    def total_attenuation_synthesis(self, lat, lon, f, el, p, D, Ns, Ts=1,
+                                    tau=45, hs=None, eta=0.65, rho=None,
+                                    H=None, P=None, hL=1000,
+                                    return_contributions=False):
+        return self.instance.total_attenuation_synthesis(
+            lat, lon, f, el, p, D, Ns, Ts=Ts, tau=tau, hs=hs, eta=eta, rho=rho,
+            H=H, P=P, hL=hL, return_contributions=return_contributions)
+
+    def scintillation_attenuation_synthesis(self, Ns, f_c=0.1, Ts=1):
+        return self.instance.scintillation_attenuation_synthesis(
+                Ns, f_c=f_c, Ts=Ts)
+
+    def cloud_liquid_water_synthesis(self, lat, lon, Ns, Ts=1, n=None):
+        return self.instance.cloud_liquid_water_synthesis(
+                lat, lon, Ns, Ts=Ts, n=n)
+
+    def integrated_water_vapour_synthesis(self, lat, lon, Ns, Ts=1, n=None):
+        return self.instance.integrated_water_vapour_synthesis(
+                lat, lon, Ns, Ts=Ts, n=n)
 
 
 class _ITU1853_1():
@@ -62,27 +84,51 @@ def __init__(self):
         self.month = 2
         self.link = 'https://www.itu.int/rec/R-REC-P.1853-1-201202-I/en'
 
-    def rain_attenuation_synthesis(self, lat, lon, f, el, hs, T, Ts=1, n=None):
+    @classmethod
+    def rain_attenuation_synthesis(self, lat, lon, f, el, hs, Ns, tau=45,
+                                   Ts=1, n=None):
         """
         For Earth-space paths, the time series synthesis method is valid for
         frequencies between 4 GHz and 55 GHz and elevation angles between
         5 deg and 90 deg.
         """
 
-        P_rain = rain_percentage_probability(lat, lon).\
+        # Step A1: Determine Prain (% of time), the probability of rain on the
+        # path. Prain can be well approximated as P0(lat, lon)
+        P_rain = rainfall_probability(lat, lon).\
             to(u.dimensionless_unscaled).value
 
-        sigma, m = fit_rain_attenuation_to_lognormal(
-            lat, lon, f, el, hs, P_rain * 100)
+        # Step A2: Construct the set of pairs [Pi, Ai] where Pi (% of time) is
+        # the probability the attenuation Ai (dB) is exceeded where Pi < P_K
+        p_i = np.array([0.01, 0.02, 0.03, 0.05,
+                        0.1, 0.2, 0.3, 0.5, 1, 2, 3, 5, 10])
+        Pi = np.array([p for p in p_i if p < P_rain * 100], dtype=np.float)
+        Ai = np.array([0 for p in p_i if p < P_rain * 100], dtype=np.float)
+
+        for i, p in enumerate(Pi):
+            Ai[i] = rain_attenuation(lat, lon, f, el, hs, p, tau=tau).value
+
+        # Step A3: Transform the set of pairs [Pi, Ai] to [Q^{-1}(Pi/P_k),
+        # ln(Ai)]
+        Q = stats.norm.ppf((Pi / 100))
+        lnA = np.log(Ai)
+
+        # Step A4: Determine the variables sigma_lna, m_lna by performing a
+        # least-squares fit to lnAi = sigma_lna Q^{-1}(Pi/P_k) + m_lna
+        m, sigma = np.linalg.lstsq(np.vstack([np.ones(len(Q)), Q]).T,
+                                   lnA, rcond=None)[0]
 
         # Step B: Set the low-pass filter parameter
-        beta = 2e-5
+        beta = 2e-4
         # Step C: compute the attenuation offset
         A_offset = np.exp(m + sigma * stats.norm.ppf(P_rain))
         # Step D: Time series synthesis
         # D1: Synthesize a white Gaussian noise time series
         if n is None:
-            n = np.random.normal(0, 1, (T * Ts + 2e5))
+            n = np.random.normal(0, 1, int(Ns * Ts + 2e5))[::Ts]
+            discard_samples = True
+        else:
+            discard_samples = False
 
         # D2, D3 : Filter the noise time series with a recursive low-pass
         # filter
@@ -94,11 +140,12 @@ def rain_attenuation_synthesis(self, lat, lon, f, el, hs, T, Ts=1, n=None):
         A_rain = np.maximum(Y_rain - A_offset, 0)
 
         # D6: Discard the first 200 000 samples from the synthesized
-        A_rain = A_rain[200000::Ts]
+        if discard_samples:
+            A_rain = A_rain[200000:]
 
-        return A_rain
+        return A_rain.flatten()
 
-    def fftnoise(f):
+    def fftnoise(self, f):
         f = np.array(f, dtype='complex')
         Np = (len(f) - 1) // 2
         phases = np.random.rand(Np) * 2 * np.pi
@@ -107,20 +154,21 @@ def fftnoise(f):
         f[-1:-1 - Np:-1] = np.conj(f[1:Np + 1])
         return np.fft.ifft(f).real
 
-    def scintillation_synthesis(self, T, f_c=0.1, Ts=1, f_c2=0.01):
+    @classmethod
+    def scintillation_attenuation_synthesis(self, Ns, f_c=0.1, Ts=1):
         """
         For Earth-space paths, the time series synthesis method is valid for
         frequencies between 4 GHz and 55 GHz and elevation angles between
         5 deg and 90 deg.
         """
-        freqs = np.abs(np.fft.fftfreq(2 * int(T + 2e5), 1 / Ts))
+        freqs = np.abs(np.fft.fftfreq(2 * int(Ns + 2e5), 1 / Ts))
         H_f = np.where(freqs <= f_c, 1, 10 **
                        ((np.log10(freqs) - np.log10(f_c)) * (-8 / 3)))
-        H_f = H_f[0:int(T + 2e5)]
+        H_f = H_f[0:int(Ns + 2e5)]
         sci = self.fftnoise(np.fft.fftshift(H_f))
-        return sci
+        return sci[200000:].flatten()
 
-    def cloud_liquid_water_synthesis(self, lat, lon, T, Ts=1, n=None):
+    def cloud_liquid_water_synthesis(self, lat, lon, Ns, Ts=1, n=None):
         """
         """
 
@@ -128,7 +176,7 @@ def cloud_liquid_water_synthesis(self, lat, lon, T, Ts=1, n=None):
         m, sigma, Pcwl = lognormal_approximation_coefficient(lat, lon)
         m = m.value
         sigma = sigma.value
-        Pcwl = Pcwl.value
+        Pcwl = Pcwl.value / 100
 
         # Step B: Low pass filter parameters
         beta_1 = 7.17e-4
@@ -142,7 +190,11 @@ def cloud_liquid_water_synthesis(self, lat, lon, T, Ts=1, n=None):
         # Step D: Time series synthesis
         # Step D1: Synthesize a white Gaussian noise time series
         if n is None:
-            n = np.random.rand((T + 5e5) // Ts)
+            n = np.random.normal(0, 1, (Ns * Ts + 5e5))[::Ts]
+            discard_samples = True
+        else:
+            discard_samples = False
+
         # Step D3: Filter the noise time series, with two recursive low-pass
         # filters
         rho_1 = np.exp(-beta_1 * Ts)
@@ -156,23 +208,29 @@ def cloud_liquid_water_synthesis(self, lat, lon, T, Ts=1, n=None):
             1 - 1 / Pcwl * stats.norm.sf(G_c))), 0)
 
         # D6: Discard the first 500 000 samples from the synthesized
-        L = L[5e5:]
+        if discard_samples:
+            L = L[500000:]
 
-        return L
+        return L.flatten()
 
-    def integrated_water_vapour_synthesis(self, lat, lon, T, Ts=1, n=None):
+    def integrated_water_vapour_coefficients(self, lat, lon):
         # A Estimation of κ and λ
         ps = np.array([0.1, 0.2, 0.3, 0.5, 1, 2, 3, 5, 10, 20, 30, 50])
-        Vi = np.array([total_water_vapour_content(lat, lon, p_i)
+        Vi = np.array([total_water_vapour_content(lat, lon, p_i).value
                        for p_i in ps])
 
-        ln_lnPi = np.log(- np.log(ps))
+        ln_lnPi = np.log(- np.log(ps / 100))
         ln_Vi = np.log(Vi)
 
         a, b = np.linalg.lstsq(np.vstack([ln_Vi, np.ones(len(ln_Vi))]).T,
                                ln_lnPi)[0]
         kappa = a
         lambd = np.exp(-b / a)
+        return kappa, lambd
+
+    def integrated_water_vapour_synthesis(self, lat, lon, Ns, Ts=1, n=None):
+        # A Estimation of κ and λ
+        kappa, lambd = self.integrated_water_vapour_coefficients(lat, lon)
 
         # B Low-pass filter parameter
         beta_V = 3.24e-6
@@ -180,7 +238,11 @@ def integrated_water_vapour_synthesis(self, lat, lon, T, Ts=1, n=None):
         # Step C: Time series synthesis
         # Step C1: Synthesize a white Gaussian noise time series
         if n is None:
-            n = np.random.rand((T + 5e6) // Ts)
+            n = np.random.normal(0, 1, (Ns * Ts + 5000000))[::Ts]
+            discard_samples = True
+        else:
+            discard_samples = False
+
         # Step C3: Filter the noise time series, with two recursive low-pass
         # filters
         rho = np.exp(-beta_V * Ts)
@@ -188,11 +250,16 @@ def integrated_water_vapour_synthesis(self, lat, lon, T, Ts=1, n=None):
         # Step C4: Compute Compute V(kTs),
         V = lambd * (- np.log10(stats.norm.sf(G_v)))**(1 / kappa)
         # Step C5: Discard the first 5 000 000 samples from the synthesized
-        V = V[5e6:]
+        if discard_samples:
+            V = V[5000000:]
 
-        return V
+        return V.flatten()
 
-    def total_attenuation_synthesis(self, lat, lon, f, el, p, hs, T, Ts=1):
+    def total_attenuation_synthesis(self, lat, lon, f, el, p, D, Ns, Ts=1,
+                                    hs=None, tau=45, eta=0.65, rho=None,
+                                    H=None, P=None, hL=1000,
+                                    return_contributions=False):
+        t_disc = int(5e6)
         # Step A Correlation coefficients:
         C_RC = 1
         C_CV = 0.8
@@ -207,50 +274,68 @@ def a_Enhanc(p):
                 np.log10(p)**2 - 1.258 * np.log10(p) + 2.672
 
         # Step C1-C3:
-        n_R = np.random.rand((T + 5e6) // Ts)
-        n_L0 = np.random.rand((T + 5e6) // Ts)
-        n_V0 = np.random.rand((T + 5e6) // Ts)
+        n_R = np.random.normal(0, 1, int((Ns * Ts + t_disc)))
+        n_L0 = np.random.normal(0, 1, int((Ns * Ts + t_disc)))
+        n_V0 = np.random.normal(0, 1, int((Ns * Ts + t_disc)))
+
         # Step C4-C5:
         n_L = C_RC * n_R + np.sqrt(1 - C_RC**2) * n_L0
         n_V = C_CV * n_L + np.sqrt(1 - C_CV**2) * n_V0
 
         # Step C6: Compute the rain attenuation time series
+        if hs is None:
+            hs = topographic_altitude(lat, lon)
         Ar = self.rain_attenuation_synthesis(
-            lat, lon, f, el, hs, T + 5e6, n=n_R)
-        Ar = Ar[5e6:]
+            lat, lon, f, el, hs, Ns, Ts=1, tau=tau, n=n_R)
+        Ar = Ar[t_disc:]
 
         # Step C7: Compute the cloud integrated liquid water content time
         # series
-        L = self.cloud_liquid_water_synthesis(lat, lon, T + 5e6, n=n_L)
+        L = self.cloud_liquid_water_synthesis(lat, lon, Ns, Ts=1, n=n_L)
+        L = L[t_disc:]
         Ac = L * \
             specific_attenuation_coefficients(f, T=0) / np.sin(np.deg2rad(el))
+        Ac = Ac.flatten()
 
         # Step C9: Identify time stamps where A_R > 0 L > 1
-        idx = np.where(Ar > 0 & L > 1)[0]
-        idx_no = np.where(Ar <= 0 | L <= 1)[0]
+        idx = np.where(np.logical_and(Ar > 0, L > 1))[0]
+        idx_no = np.where(np.logical_not(
+                    np.logical_and(Ar > 0, L > 1)))[0]
+
         # Step C10: Discard the previous values of Ac and re-compute them by
         # linear interpolation vs. time starting from the non-discarded cloud
         # attenuations values
         Ac[idx] = np.interp(idx, idx_no, Ac[idx_no])
 
         # Step C11: Compute the integrated water vapour content time series
-        V = self.integrated_water_vapour_synthesis(lat, lon, T, n=n_V)
+        V = self.integrated_water_vapour_synthesis(lat, lon, Ns, Ts=1, n=n_V)
+        V = V[t_disc:]
 
         # Step C12: Convert the integrated water vapour content time series
         # V into water vapour attenuation time series AV(kTs)
-        Av = zenit_water_vapour_attenuation(lat, lon, p, f, V_t=V)
+        Av = zenit_water_vapour_attenuation(lat, lon, p, f, V_t=V).value
 
         # Step C13: Compute the mean annual temperature Tm for the location of
         # interest using experimental values if available.
-        Tm = surface_mean_temperature(lat, lon)
+        Tm = surface_mean_temperature(lat, lon).value
 
         # Step C14: Convert the mean annual temperature Tm into mean annual
         # oxygen attenuation AO following the method recommended in
         # Recommendation ITU-R P.676.
-        Ao = 1
+        if P is None:
+            P = standard_pressure(hs).value
+
+        if rho is None:
+            rho = standard_water_vapour_density(hs).value
+
+        e = Tm * rho / 216.7
+        go = gamma0_exact(f, P, rho, Tm).value
+        ho, hw = slant_inclined_path_equivalent_height(f, P + e).value
+        Ao = ho * go * np.ones_like(Ar)
 
         # Step C15: Synthesize unit variance scintillation time series
-        sci_0 = self.scintillation_synthesis(T)
+        sci_0 = self.scintillation_attenuation_synthesis(Ns * Ts, Ts=1)
+
         # Step C16: Compute the correction coefficient time series Cx(kTs) in
         # order to distinguish between scintillation fades and enhancements:
         Q_sci = 100 * stats.norm.sf(sci_0)
@@ -258,21 +343,56 @@ def a_Enhanc(p):
 
         # Step C17: Transform the integrated water vapour content time series
         # V(kTs) into the Gamma distributed time series Z(kTs) as follows:
-        Z = scipy.stats.gamma.ppf(1 - np.exp(-(V / lambd)**kappa), 10, 0.1)
+        kappa, lambd = self.integrated_water_vapour_coefficients(lat, lon)
+        Z = stats.gamma.ppf(1 - np.exp(-(V / lambd)**kappa), 10, 0.1)
 
         # Step C18: Compute the scintillation standard deviation σ following
         # the method recommended in Recommendation ITU-R P.618.
-        sigma = 2
+        sigma = scintillation_attenuation_sigma(lat, lon, f, el, p, D, eta, Tm,
+                                                H, P, hL).value
 
         # Step C19: Compute the scintillation time series sci:
-        sci = np.where(Ar > 1, sigma * sci_0 * C_x * Z * Ar ** (5 / 12),
-                       sigma * sci_0 * C_x * Z)
+        As = np.where(Ar > 1, sigma * sci_0 * C_x * Z * Ar ** (5 / 12),
+                      sigma * sci_0 * C_x * Z)
 
         # Step C20: Compute total tropospheric attenuation time series A(kTs)
         # as follows:
-        A = Ar + Ac + Av + Ao + sci
+        A = Ar + Ac + Av + Ao + As
+
+        if return_contributions:
+            return (Ao + Av)[::Ts], Ac[::Ts], Ar[::Ts], As[::Ts], A[::Ts]
+        else:
+            return A[::Ts]
+
+
+class _ITU1853_0():
+
+    def __init__(self):
+        self.__version__ = 0
+        self.year = 2009
+        self.month = 10
+        self.link = 'https://www.itu.int/rec/R-REC-P.1853-0-200910-I/en'
+
+    def rain_attenuation_synthesis(*args, **kwargs):
+        return _ITU1853_1.rain_attenuation_synthesis(*args, **kwargs)
+
+    def scintillation_attenuation_synthesis(self, *args, **kwargs):
+        return _ITU1853_1.scintillation_attenuation_synthesis(*args, **kwargs)
+
+    def cloud_liquid_water_synthesis(*args, **kwargs):
+        raise NotImplementedError(
+            "Recommendation ITU-R P.1853 does not specify a method to compute "
+            "time series for the cloud liquid water content.")
+
+    def integrated_water_vapour_synthesis(*args, **kwargs):
+        raise NotImplementedError(
+            "Recommendation ITU-R P.1853 does not specify a method to compute "
+            "time series for the water vapour content.")
 
-        return A
+    def total_attenuation_synthesis(*args, **kwargs):
+        raise NotImplementedError(
+            "Recommendation ITU-R P.1853 does not specify a method to compute "
+            "time series for the total atmospheric attenuation.")
 
 
 __model = __ITU1853()
@@ -288,12 +408,7 @@ def change_version(new_version):
     new_version : int
         Number of the version to use.
         Valid values are:
-        * P.1853-1 (08/94) (Superseded)
-        * P.1853-2 (10/99) (Superseded)
-        * P.1853-3 (02/01) (Superseded)
-        * P.1853-4 (04/03) (Superseded)
-        * P.1853-5 (08/07) (Superseded)
-        * P.1853-6 (02/12) (Current version)
+        * P.1853-1 (02/12) (Current version)
     """
     global __model
     __model = __ITU1853(new_version)
@@ -307,9 +422,9 @@ def get_version():
     return __model.__version__
 
 
-def rain_attenuation_synthesis(lat, lon, f, el, hs, T, Ts=1, n=None):
+def rain_attenuation_synthesis(lat, lon, f, el, hs, Ns, Ts=1, tau=45, n=None):
     """
-    A method to compute the rainfall rate exceeded for p% of the average year
+    A method to generate a synthetic time series of rain attenuation values.
 
 
     Parameters
@@ -327,10 +442,13 @@ def rain_attenuation_synthesis(lat, lon, f, el, hs, T, Ts=1, n=None):
         the earth station height above mean sea level is not available, an
         estimate is obtained from the maps of topographic altitude
         given in Recommendation ITU-R P.1511.
-    T : int
+    Ns : int
         Number of samples
     Ts : int
         Time step between consecutive samples (seconds)
+    tau : number, optional
+        Polarization tilt angle relative to the horizontal (degrees)
+        (tau = 45 deg for circular polarization). Default value is 45
     n : list, np.array, optional
         Additive White Gaussian Noise used as input for the
 
@@ -354,15 +472,47 @@ def rain_attenuation_synthesis(lat, lon, f, el, hs, T, Ts=1, n=None):
     hs = prepare_quantity(
         hs, u.km, 'Heigh above mean sea level of the earth station')
     Ts = prepare_quantity(f, u.second, 'Time step between samples')
-    val = __model.rain_attenuation_synthesis(lat, lon, f, el, hs, T, Ts, n)
+    val = __model.rain_attenuation_synthesis(lat, lon, f, el, hs, Ns,
+                                             T=Ts, tau=tau, n=n)
     return val * u.dB
 
 
-def unavailability_from_rainfall_rate(lat, lon, R):
+def scintillation_attenuation_synthesis(self, Ns, f_c=0.1, Ts=1):
     """
-    A method to estimate the percentage of time of the average year that a given
-    rainfall rate (R) is exceeded. This method calls successively to the
-    `rainfall_rate` method and interpolates its value.
+    A method to generate a synthetic time series of scintillation attenuation
+    values.
+
+
+    Parameters
+    ----------
+    Ns : int
+        Number of samples
+    f_c : float
+        Cut-off frequency for the low pass filter
+    Ts : int
+        Time step between consecutive samples (seconds)
+
+    Returns
+    -------
+    sci_att: numpy.ndarray
+        Synthesized scintilation attenuation time series (dB)
+
+
+    References
+    ----------
+    [1] Characteristics of precipitation for propagation modelling
+    https://www.itu.int/rec/R-REC-P.1853/en
+    """
+    global __model
+
+    val = __model.scintillation_attenuation_synthesis(Ns, f_c, Ts)
+    return val * u.dB
+
+
+def integrated_water_vapour_synthesis(self, lat, lon, Ns, Ts=1, n=None):
+    """ The time series synthesis method generates a time series that
+    reproduces the spectral characteristics and the distribution of water
+    vapour content.
 
 
     Parameters
@@ -371,14 +521,53 @@ def unavailability_from_rainfall_rate(lat, lon, R):
         Latitudes of the receiver points
     lon : number, sequence, or numpy.ndarray
         Longitudes of the receiver points
-    R : number, sequence, or numpy.ndarray
-        Rainfall rate (mm/h)
+    Ns : int
+        Number of samples
+    Ts : int
+        Time step between consecutive samples (seconds)
+    n : list, np.array, optional
+        Additive White Gaussian Noise used as input for the
+
+    Returns
+    -------
+    L: numpy.ndarray
+        Synthesized water vapour content time series (kg/m2)
+
 
+    References
+    ----------
+    [1] Characteristics of precipitation for propagation modelling
+    https://www.itu.int/rec/R-REC-P.1853/en
+    """
+    global __model
+
+    lon = np.mod(lon, 360)
+    val = __model.integrated_water_vapour_synthesis(lat, lon, Ns, Ts, n)
+    return val * u.kg / u.m**2
+
+
+def cloud_liquid_water_synthesis(self, lat, lon, Ns, Ts=1, n=None):
+    """ The time series synthesis method generates a time series that
+    reproduces the spectral characteristics, rate of change and duration
+    statistics of cloud liquid content events.
+
+    Parameters
+    ----------
+    lat : number, sequence, or numpy.ndarray
+        Latitudes of the receiver points
+    lon : number, sequence, or numpy.ndarray
+        Longitudes of the receiver points
+    Ns : int
+        Number of samples
+    Ts : int
+        Time step between consecutive samples (seconds)
+    n : list, np.array, optional
+        Additive White Gaussian Noise used as input for the
 
     Returns
     -------
-    p: numpy.ndarray
-        Rainfall rate exceeded for p% of the average year
+    V: numpy.ndarray
+        Synthesized cloud liquid water time series (mm)
 
 
     References
@@ -388,4 +577,98 @@ def unavailability_from_rainfall_rate(lat, lon, R):
     """
     global __model
 
-    # TODO: write this function
+    lon = np.mod(lon, 360)
+    val = __model.cloud_liquid_water_synthesis(lat, lon, Ns, Ts, n)
+    return val * u.mm
+
+
+def total_attenuation_synthesis(lat, lon, f, el, p, D, Ns, Ts=1, hs=None,
+                                tau=45, eta=0.65,  rho=None, H=None, P=None,
+                                hL=1000, return_contributions=False):
+    """ The time series synthesis method generates a time series that
+    reproduces the spectral characteristics, rate of change and duration
+    statistics of the total atmospheric attenuation events.
+
+    The time series is obtained considering the contributions of gaseous,
+    cloud, rain, and scintillation attenuation.
+
+    Parameters
+    ----------
+    lat : number
+        Latitudes of the receiver points
+    lon : number
+        Longitudes of the receiver points
+    f : number or Quantity
+        Frequency (GHz)
+    el : number
+        Elevation angle (degrees)
+    p : number
+        Percetage of the time the rain attenuation value is exceeded.
+    D: number or Quantity
+        Physical diameter of the earth-station antenna (m)
+    Ns : int
+        Number of samples
+    Ts : int
+        Time step between consecutive samples (seconds)
+    tau : number, optional
+        Polarization tilt angle relative to the horizontal (degrees)
+        (tau = 45 deg for circular polarization). Default value is 45
+    hs : number, sequence, or numpy.ndarray, optional
+        Heigh above mean sea level of the earth station (km). If local data for
+        the earth station height above mean sea level is not available, an
+        estimate is obtained from the maps of topographic altitude
+        given in Recommendation ITU-R P.1511.
+    eta: number, optional
+        Antenna efficiency. Default value 0.5 (conservative estimate)
+    rho : number or Quantity, optional
+        Water vapor density (g/m3). If not provided, an estimate is obtained
+        from Recommendation Recommendation ITU-R P.836.
+    H: number, sequence, or numpy.ndarray, optional
+        Average surface relative humidity (%) at the site. If None, uses the
+        ITU-R P.453 to estimate the wet term of the radio refractivity.
+    P: number, sequence, or numpy.ndarray, optional
+        Average surface pressure (hPa) at the site. If None, uses the
+        ITU-R P.453 to estimate the wet term of the radio refractivity.
+    hL : number, optional
+        Height of the turbulent layer (m). Default value 1000 m
+    return_contributions: bool, optional
+        Determines whether individual contributions from gases, rain, clouds
+        and scintillation are returned in addition ot the total attenuation
+        (True), or just the total atmospheric attenuation (False).
+        Default is False
+
+    Returns
+    ---------
+    A : Quantity
+        Synthesized total atmospheric attenuation time series (dB)
+
+    Ag, Ac, Ar, As, A : tuple
+        Synthesized Gaseous, Cloud, Rain, Scintillation contributions to total
+        attenuation time series, and synthesized total attenuation time seires
+        (dB).
+
+    References
+    ----------
+    [1] Characteristics of precipitation for propagation modelling
+    https://www.itu.int/rec/R-REC-P.1853/en
+    """
+    global __model
+
+    f = prepare_quantity(f, u.GHz, 'Frequency')
+    el = prepare_quantity(el, u.deg, 'Elevation angle')
+    D = prepare_quantity(D, u.m, 'Antenna diameter')
+    hs = prepare_quantity(
+        hs, u.km, 'Heigh above mean sea level of the earth station')
+    eta = prepare_quantity(eta, u.one, 'Antenna efficiency')
+    rho = prepare_quantity(rho, u.g / u.m**3, 'Water vapor density')
+    H = prepare_quantity(H, u.percent, 'Average surface relative humidity')
+    P = prepare_quantity(P, u.hPa, 'Average surface pressure')
+    hL = prepare_quantity(hL, u.m, 'Height of the turbulent layer')
+
+    val = __model.total_attenuation_synthesis(
+            lat, lon, f, el, p, D, Ns, Ts=Ts, tau=tau, hs=hs, eta=eta, rho=rho,
+            H=H, P=P, hL=hL, return_contributions=return_contributions)
+    if return_contributions:
+        return tuple([v * u.dB for v in val])
+    else:
+        return val * u.dB
diff --git a/itur/models/itu530.py b/itur/models/itu530.py
index 87b5bf9..3a7ae9d 100644
--- a/itur/models/itu530.py
+++ b/itur/models/itu530.py
@@ -6,6 +6,7 @@
 import numpy as np
 import os
 from astropy import units as u
+from scipy.optimize import bisect
 
 from itur.models.itu453 import DN65
 from itur.models.itu837 import rainfall_rate
@@ -36,8 +37,10 @@ class __ITU530():
     # This is an abstract class that contains an instance to a version of the
     # ITU-R P.530 recommendation.
 
-    def __init__(self, version=16):
-        if version == 16:
+    def __init__(self, version=17):
+        if version == 17:
+            self.instance = _ITU530_17()
+        elif version == 16:
             self.instance = _ITU530_16()
         else:
             raise ValueError(
@@ -84,16 +87,17 @@ def XPD_outage_precipitation(self, lat, lon, d, f, el, C0_I, tau=45,
                                                       tau, U0, XPIF)
 
 
-class _ITU530_16():
+class _ITU530_17():
 
     def __init__(self):
-        self.__version__ = 16
-        self.year = 2015
-        self.month = 7
-        self.link = 'https://www.itu.int/rec/R-REC-P.530-16-201507-S/en'
+        self.__version__ = 17
+        self.year = 2017
+        self.month = 12
+        self.link = 'https://www.itu.int/rec/R-REC-P.530-17-201712-S/en'
 
         self._s_a = {}
 
+    @classmethod
     def s_a(self, lat, lon):
         """ Standard deviation of terrain heights (m) within a 110 km × 110 km
         area with a 30 s resolution (e.g. the Globe “gtopo30” data).
@@ -111,8 +115,9 @@ def s_a(self, lat, lon):
             np.array([lat.ravel(), lon.ravel()]).T).reshape(lat.shape)
 
     ###########################################################################
-    ###                             Section 2.2                             ###
+    #                               Section 2.2                               #
     ###########################################################################
+    @classmethod
     def fresnel_ellipse_radius(self, d1, d2, f):
         """ Implementation of 'fresnel_ellipse_radius' method for recommendation
         ITU-P R.530-16. See documentation for function
@@ -120,6 +125,7 @@ def fresnel_ellipse_radius(self, d1, d2, f):
         """
         return 17.3 * np.sqrt(d1 * d2 / (f * (d1 + d2)))
 
+    @classmethod
     def diffraction_loss(self, d1, d2, h, f):
         """ Implementation of 'diffraction_loss' method for recommendation
         ITU-P R.530-16. See documentation for function
@@ -130,8 +136,9 @@ def diffraction_loss(self, d1, d2, h, f):
         return Ad
 
     ###########################################################################
-    ###                             Section 2.3                             ###
+    #                               Section 2.3                               #
     ###########################################################################
+    @classmethod
     def multipath_loss_for_A(self, lat, lon, h_e, h_r, d, f, A):
         """ Implementation of 'multipath_loss_for_A' method for recommendation
         ITU-P R.530-16. See documentation for function
@@ -150,15 +157,16 @@ def multipath_loss_for_A(self, lat, lon, h_e, h_r, d, f, A):
         # Eq. 5 [mrad]
         e_p = np.abs(h_r - h_e) / d
 
-        # Step 3: For detailed link design applications calculate the percentage
-        # of time (p_W) that fade depth A (dB) is exceeded in the average worst
-        # month
+        # Step 3: For detailed link design applications calculate the
+        # percentage of time (p_W) that fade depth A (dB) is exceeded in the
+        # average worst month
         h_L = np.minimum(h_e, h_r)
         p_W = K * d**3.4 (1 + e_p)**-1.03 * f**0.8 * \
             10**(-0.00076 * h_L - A / 10)
         # Eq. 7 [%]
         return p_W
 
+    @classmethod
     def multipath_loss(self, lat, lon, h_e, h_r, d, f, A):
         """ Implementation of 'multipath_loss' method for recommendation
         ITU-P R.530-16. See documentation for function
@@ -169,9 +177,9 @@ def multipath_loss(self, lat, lon, h_e, h_r, d, f, A):
         p0 = self.multipath_loss_for_A(
             lat, lon, h_e, h_r, d, f, 0)   # Eq. 10 [%]
 
-        # Step 2: Calculate the value of fade depth, At, at which the transition
-        # occurs between the deep-fading distribution and the shallow-fading
-        # distribution
+        # Step 2: Calculate the value of fade depth, At, at which the
+        # transition occurs between the deep-fading distribution and the
+        # shallow-fading distribution
         At = 25 + 1.2 * np.log10(p0)                        # Eq. 12 [dB]
 
         # Step 3: Calculate the percentage of time that A is exceeded in the
@@ -179,8 +187,9 @@ def multipath_loss(self, lat, lon, h_e, h_r, d, f, A):
         def step_3b(p_0, At, A):
             p_t = p_0 * 10 ** (-At / 10)
             qa_p = -20 * np.log10(-np.log((100 - p_t) / 100)) / At
-            q_t = (qa_p - 2) / (1 + 0.3 * 10 ** (-At / 20) * 10 **
-                                (-0.016 * At)) - 4.3 * (10**(-At / 20) + At / 800)
+            q_t = ((qa_p - 2) /
+                   (1 + 0.3 * 10 ** (-At / 20) * 10 ** (-0.016 * At)) -
+                   4.3 * (10**(-At / 20) + At / 800))
             q_a = 2 + (1 + 0.3 * 10**(-A / 20)) * (10**(-0.016 * A)) *\
                 (q_t + 4.3 * (10**(-A / 20 + A / 800)))
             p_W = 100 * (1 - np.exp(-10 ** (-q_a * A / 20)))
@@ -191,8 +200,9 @@ def step_3b(p_0, At, A):
         return p_W
 
     ###########################################################################
-    ###                             Section 2.4                             ###
+    #                               Section 2.4                               #
     ###########################################################################
+    @classmethod
     def rain_attenuation(self, lat, lon, d, f, el, p, tau=45, R001=None):
         """ Implementation of 'rain_attenuation' method for recommendation
         ITU-P R.530-16. See documentation for function
@@ -207,7 +217,7 @@ def rain_attenuation(self, lat, lon, d, f, el, p, tau=45, R001=None):
         # frequency, polarization and rain rate of interest using
         # Recommendation ITU-R P.838
         gammar = rain_specific_attenuation(R001, f, el, tau).value
-        _, alpha = rain_specific_attenuation_coefficients(f, el, tau).value
+        _, alpha = rain_specific_attenuation_coefficients(f, el, tau)
 
         # Step 3: Compute the effective path length, deff, of the link by
         # multiplying the actual path length d by a distance factor r
@@ -232,8 +242,8 @@ def rain_attenuation(self, lat, lon, d, f, el, p, tau=45, R001=None):
 
     def inverse_rain_attenuation(
             self, lat, lon, d, f, el, Ap, tau=45, R001=None):
-        """ Implementation of 'inverse_rain_attenuation' method for recommendation
-        ITU-P R.530-16. See documentation for function
+        """ Implementation of 'inverse_rain_attenuation' method for
+        recommendation ITU-P R.530-16. See documentation for function
         'ITUR530.inverse_rain_attenuation'
         """
         # Step 1: Obtain the rain rate R0.01 exceeded for 0.01% of the time
@@ -264,9 +274,12 @@ def inverse_rain_attenuation(
         C2 = 0.855 * C0 + 0.546 * (1 - C0)
         C3 = 0.139 * C0 + 0.043 * (1 - C0)
 
-        Ap = A001 * C1 * p ** (- (C2 + C3 * np.log10(p)))
-        return Ap
+        def func_bisect(p):
+            return A001 * C1 * p ** (- (C2 + C3 * np.log10(p))) - Ap
 
+        return bisect(func_bisect, 0, 100)
+
+    @classmethod
     def rain_event_count(self, lat, lon, d, f, el, A, tau=45, R001=None):
         """ Implementation of 'rain_event_count' method for recommendation
         ITU-P R.530-16. See documentation for function
@@ -282,9 +295,9 @@ def rain_event_count(self, lat, lon, d, f, el, A, tau=45, R001=None):
         return N10s
 
     ###########################################################################
-    ###                              Section 4                              ###
+    #                                Section 4                                #
     ###########################################################################
-
+    @classmethod
     def XPD_outage_clear_air(self, lat, lon, h_e, h_r,
                              d, f, XPD_g, C0_I, XPIF=0):
         """ Implementation of 'XPD_outage_clear_air' method for recommendation
@@ -311,6 +324,7 @@ def XPD_outage_clear_air(self, lat, lon, h_e, h_r,
         P_XP = P0 * 10 ** (- M_XPD / 10)                       # Eq. 106 [%]
         return P_XP
 
+    @classmethod
     def XPD_outage_precipitation(self, lat, lon, d, f, el, C0_I, tau=45,
                                  U0=15, XPIF=0):
         """ Implementation of 'XPD_outage_precipitation' method for recommendation
@@ -335,6 +349,59 @@ def XPD_outage_precipitation(self, lat, lon, d, f, el, C0_I, tau=45,
         return P_XPR
 
 
+class _ITU530_16():
+
+    def __init__(self):
+        self.__version__ = 16
+        self.year = 2015
+        self.month = 7
+        self.link = 'https://www.itu.int/rec/R-REC-P.530-16-201507-S/en'
+
+        self._s_a = {}
+
+    def s_a(self, *args, **kwargs):
+        return _ITU530_17.s_a(*args, **kwargs)
+
+    ###########################################################################
+    #                               Section 2.2                               #
+    ###########################################################################
+    def fresnel_ellipse_radius(self, *args, **kwargs):
+        return _ITU530_17.fresnel_ellipse_radius(*args, **kwargs)
+
+    def diffraction_loss(self, *args, **kwargs):
+        return _ITU530_17.diffraction_loss(*args, **kwargs)
+
+    ###########################################################################
+    #                               Section 2.3                               #
+    ###########################################################################
+    def multipath_loss_for_A(self, *args, **kwargs):
+        return _ITU530_17.multipath_loss_for_A(*args, **kwargs)
+
+    def multipath_loss(*args, **kwargs):
+        return _ITU530_17.multipath_loss(*args, **kwargs)
+
+    ###########################################################################
+    #                               Section 2.4                               #
+    ###########################################################################
+    def rain_attenuation(self, *args, **kwargs):
+        return _ITU530_17.rain_attenuation(*args, **kwargs)
+
+    def inverse_rain_attenuation(self, *args, **kwargs):
+        return _ITU530_17.inverse_rain_attenuation(*args, **kwargs)
+
+    def rain_event_count(self, *args, **kwargs):
+        return _ITU530_17.rain_event_count(*args, **kwargs)
+
+    ###########################################################################
+    #                                Section 4                                #
+    ###########################################################################
+    def XPD_outage_clear_air(self, *args, **kwargs):
+        return _ITU530_17.XPD_outage_clear_air(*args, **kwargs)
+
+    def XPD_outage_precipitation(self, *args, **kwargs):
+        return _ITU530_17.XPD_outage_precipitation(*args, **kwargs)
+
+
 __model = __ITU530()
 
 
@@ -564,10 +631,10 @@ def multipath_loss(lat, lon, h_e, h_r, d, f, A):
 
 def rain_attenuation(lat, lon, d, f, el, p, tau=45, R001=None):
     """ Estimate long-term statistics of rain attenuation. Attenuation can also
-    occur as a result of absorption and scattering by such hydrometeors as rain,
-    snow, hail and fog. Although rain attenuation can be ignored at frequencies
-    below about 5 GHz, it must be included in design calculations at higher
-    frequencies, where its importance increases rapidly.
+    occur as a result of absorption and scattering by such hydrometeors as
+    rain, snow, hail and fog. Although rain attenuation can be ignored at
+    frequencies below about 5 GHz, it must be included in design calculations
+    at higher frequencies, where its importance increases rapidly.
 
 
     Parameters
@@ -585,8 +652,8 @@ def rain_attenuation(lat, lon, d, f, el, p, tau=45, R001=None):
     p : number
         Percetage of the time the rain attenuation value is exceeded.
     R001: number, optional
-        Point rainfall rate for the location for 0.01% of an average year (mm/h).
-        If not provided, an estimate is obtained from Recommendation
+        Point rainfall rate for the location for 0.01% of an average year
+        (mm/h). If not provided, an estimate is obtained from Recommendation
         Recommendation ITU-R P.837. Some useful values:
             * 0.25 mm/h : Drizle
             *  2.5 mm/h : Light rain
@@ -645,8 +712,8 @@ def inverse_rain_attenuation(lat, lon, d, f, el, Ap, tau=45, R001=None):
     Ap : number
         Fade depth
     R001: number, optional
-        Point rainfall rate for the location for 0.01% of an average year (mm/h).
-        If not provided, an estimate is obtained from Recommendation
+        Point rainfall rate for the location for 0.01% of an average year
+        (mm/h). If not provided, an estimate is obtained from Recommendation
         Recommendation ITU-R P.837. Some useful values:
             * 0.25 mm/h : Drizle
             *  2.5 mm/h : Light rain
@@ -706,8 +773,8 @@ def rain_event_count(lat, lon, d, f, el, A, tau=45, R001=None):
     A : number
         Fade depth
     R001: number, optional
-        Point rainfall rate for the location for 0.01% of an average year (mm/h).
-        If not provided, an estimate is obtained from Recommendation
+        Point rainfall rate for the location for 0.01% of an average year
+        (mm/h). If not provided, an estimate is obtained from Recommendation
         Recommendation ITU-R P.837. Some useful values:
             * 0.25 mm/h : Drizle
             *  2.5 mm/h : Light rain
diff --git a/itur/models/itu618.py b/itur/models/itu618.py
index 24d20e8..18b163f 100644
--- a/itur/models/itu618.py
+++ b/itur/models/itu618.py
@@ -136,6 +136,12 @@ def scintillation_attenuation(self, lat, lon, f, el, p, D, eta,
                            excluded=[0, 1, 3, 7, 8, 9], otypes=[np.ndarray])
         return np.array(fcn(lat, lon, f, el, p, D, eta, T, H, P, hL).tolist())
 
+    def scintillation_attenuation_sigma(self, lat, lon, f, el, p, D, eta,
+                                  T, H, P, hL):
+        fcn = np.vectorize(self.instance.scintillation_attenuation_sigma,
+                           excluded=[0, 1, 3, 7, 8, 9], otypes=[np.ndarray])
+        return np.array(fcn(lat, lon, f, el, p, D, eta, T, H, P, hL).tolist())
+
     def fit_rain_attenuation_to_lognormal(self, lat, lon, f, el, hs, P_k, tau):
         fcn = np.vectorize(self.instance.fit_rain_attenuation_to_lognormal)
         return fcn(lat, lon, f, el, hs, P_k, tau)
@@ -284,7 +290,7 @@ def fit_rain_attenuation_to_lognormal(self, lat, lon, f, el, hs, P_k, tau):
 
         # Step 1: Construct the set of pairs [Pi, Ai] where Pi (% of time) is
         # the probability the attenuation Ai (dB) is exceeded where Pi < P_K
-        p_i = np.array([0.001, 0.002, 0.003, 0.005, 0.01, 0.02, 0.03, 0.05,
+        p_i = np.array([0.01, 0.02, 0.03, 0.05,
                         0.1, 0.2, 0.3, 0.5, 1, 2, 3, 5, 10])
         Pi = np.array([p for p in p_i if p < P_k], dtype=np.float)
         Ai = np.array([0 for p in p_i if p < P_k], dtype=np.float)
@@ -300,7 +306,7 @@ def fit_rain_attenuation_to_lognormal(self, lat, lon, f, el, hs, P_k, tau):
         # Step 3: Determine the variables sigma_lna, m_lna by performing a
         # least-squares fit to lnAi = sigma_lna Q^{-1}(Pi/P_k) + m_lna
         m_lna, sigma_lna = np.linalg.lstsq(np.vstack([np.ones(len(Q)), Q]).T,
-                                           lnA)[0]
+                                           lnA, rcond=None)[0]
 
         return sigma_lna, m_lna
 
@@ -432,8 +438,9 @@ def rain_cross_polarization_discrimination(self, Ap, f, el, p, tau):
         return XPD_p
 
     @classmethod
-    def scintillation_attenuation(self, lat, lon, f, el, p, D, eta=0.5, T=None,
-                                  H=None, P=None, hL=1000):
+
+    def scintillation_attenuation_sigma(self, lat, lon, f, el, p, D, eta=0.5,
+                                        T=None, H=None, P=None, hL=1000):
         # Step 1: For the value of t, calculate the saturation water vapour
         # pressure, es, (hPa), as specified in Recommendation ITU-R P.453.
         if T is not None and H is not None and P is not None:
@@ -467,7 +474,15 @@ def scintillation_attenuation(self, lat, lon, f, el, p, D, eta=0.5, T=None,
         # Step 7: Calculate the standard deviation of the signal for the
         # applicable period and propagation path:
         sigma = sigma_ref * f**(7. / 12) * g / np.sin(np.deg2rad(el))**1.2
+        return sigma
 
+    @classmethod
+    def scintillation_attenuation(self, lat, lon, f, el, p, D, eta=0.5, T=None,
+                                  H=None, P=None, hL=1000):
+        # Step 1 - 7: Calculate the standard deviation of the signal for the
+        # applicable period and propagation path:
+        sigma = self.scintillation_attenuation_sigma(lat, lon, f, el, p,
+                                                     D, eta, T, H, P, hL)
         # Step 8: Calculate the time percentage factor, a(p), for the time
         # percentage, p, in the range between 0.01% < p < 50%:
         a = -0.061 * np.log10(p)**3 + 0.072 * \
@@ -899,6 +914,77 @@ def scintillation_attenuation(lat, lon, f, el, p, D, eta=0.5, T=None,
     return prepare_output_array(val, type_output) * u.dB
 
 
+def scintillation_attenuation_sigma(lat, lon, f, el, p, D, eta=0.5, T=None,
+                                    H=None, P=None, hL=1000):
+    """
+    Calculation of the standard deviation of the amplitude of the
+    scintillations attenuation at elevation angles greater than 5° and
+    frequencies up to 20 GHz.
+
+
+    Parameters
+    ----------
+    lat : number, sequence, or numpy.ndarray
+        Latitudes of the receiver points
+    lon : number, sequence, or numpy.ndarray
+        Longitudes of the receiver points
+    f : number or Quantity
+        Frequency (GHz)
+    el : sequence, or number
+        Elevation angle (degrees)
+    p : number
+        Percetage of the time the scintillation attenuation value is exceeded.
+    D: number
+        Physical diameter of the earth-station antenna (m)
+    eta: number, optional
+        Antenna efficiency. Default value 0.5 (conservative estimate)
+    T: number, sequence, or numpy.ndarray, optional
+        Average surface ambient temperature (°C) at the site. If None, uses the
+        ITU-R P.453 to estimate the wet term of the radio refractivity.
+    H: number, sequence, or numpy.ndarray, optional
+        Average surface relative humidity (%) at the site. If None, uses the
+        ITU-R P.453 to estimate the wet term of the radio refractivity.
+    P: number, sequence, or numpy.ndarray, optional
+        Average surface pressure (hPa) at the site. If None, uses the
+        ITU-R P.453 to estimate the wet term of the radio refractivity.
+    hL : number, optional
+        Height of the turbulent layer (m). Default value 1000 m
+
+
+    Returns
+    -------
+    attenuation: Quantity
+        Attenuation due to scintillation (dB)
+
+
+    References
+    ----------
+    [1] Propagation data and prediction methods required for the design of
+    Earth-space telecommunication systems:
+    https://www.itu.int/dms_pubrec/itu-r/rec/p/R-REC-P.618-12-201507-I!!PDF-E.pdf
+    """
+    global __model
+    type_output = type(lat)
+
+    lat = prepare_input_array(lat)
+    lon = prepare_input_array(lon)
+
+    lon = np.mod(lon, 360)
+    f = prepare_quantity(f, u.GHz, 'Frequency')
+    el = prepare_quantity(prepare_input_array(el), u.deg, 'Elevation angle')
+    D = prepare_quantity(D, u.m, 'Antenna diameter')
+    eta = prepare_quantity(eta, u.one, 'Antenna efficiency')
+    T = prepare_quantity(T, u.deg_C, 'Average surface temperature')
+    H = prepare_quantity(H, u.percent, 'Average surface relative humidity')
+    P = prepare_quantity(P, u.hPa, 'Average surface pressure')
+    hL = prepare_quantity(hL, u.m, 'Height of the turbulent layer')
+
+    val = __model.scintillation_attenuation_sigma(
+        lat, lon, f, el, p, D, eta=eta, T=T, H=H, P=P, hL=hL)
+
+    return prepare_output_array(val, type_output) * u.dB
+
+
 def rain_cross_polarization_discrimination(Ap, f, el, p, tau=45):
     """
     Calculation of the cross-polarization discrimination (XPD) statistics from
diff --git a/itur/models/itu676.py b/itur/models/itu676.py
index 1bbcd7a..e73e282 100644
--- a/itur/models/itu676.py
+++ b/itur/models/itu676.py
@@ -143,12 +143,17 @@ def gaseous_attenuation_slant_path(self, f, el, rho, P, T, V_t, h, mode):
         fcn = np.vectorize(self.instance.gaseous_attenuation_slant_path)
         return fcn(f, el, rho, P, T, V_t, h, mode)
 
+    def slant_inclined_path_equivalent_height(self, f, p):
+        fcn = np.vectorize(self.instance.slant_inclined_path_equivalent_height,
+                           excluded=[0], otypes=[np.ndarray])
+        return np.array(fcn(f, p).tolist())
+
     def zenit_water_vapour_attenuation(
             self, lat, lon, p, f, V_t=None, h=None):
         # Abstract method to compute the water vapour attenuation over the
         # slant path
         fcn = np.vectorize(self.instance.zenit_water_vapour_attenuation,
-                           excluded=[0, 1], otypes=[np.ndarray])
+                           excluded=[0, 1, 4, 5], otypes=[np.ndarray])
         return np.array(fcn(lat, lon, p, f, V_t, h).tolist())
 
     def gamma_exact(self, f, p, rho, t):
@@ -322,7 +327,7 @@ def gaseous_attenuation_approximation(self, f, el, rho, P, T):
         return gamma0, gammaw
 
     @classmethod
-    def slant_inclined_path_coefficients(self, f, p):
+    def slant_inclined_path_equivalent_height(self, f, p):
         """
         """
         rp = p / 1013.25
@@ -370,7 +375,7 @@ def gaseous_attenuation_slant_path(self, f, el, rho, P, T, V_t=None,
                 f, el, rho, P, T)
 
             e = rho * T / 216.7
-            h0, hw = self.slant_inclined_path_coefficients(f, P + e)
+            h0, hw = self.slant_inclined_path_equivalent_height(f, P + e)
 
             # Use the zenit water-vapour method if the values of V_t
             # and h are provided
@@ -430,7 +435,7 @@ def gaseous_attenuation_inclined_path(
             gammaw = 0
 
         e = rho * T / 216.7
-        h0, hw = self.slant_inclined_path_coefficients(f, P + e)
+        h0, hw = self.slant_inclined_path_equivalent_height(f, P + e)
 
         if 5 < el and el < 90:
             h0_p = h0 * (np.exp(-h1 / h0) - np.exp(-h2 / h0))
@@ -662,7 +667,7 @@ def gaseous_attenuation_approximation(self, f, el, rho, P, T):
         return gamma0, gammaw
 
     @classmethod
-    def slant_inclined_path_coefficients(self, f, p):
+    def slant_inclined_path_equivalent_height(self, f, p):
         """
         """
         rp = p / 1013.0
@@ -709,7 +714,7 @@ def gaseous_attenuation_slant_path(self, f, el, rho, P, T, V_t=None,
             gamma0, gammaw = self.gaseous_attenuation_approximation(
                 f, el, rho, P, T)
             e = rho * T / 216.7
-            h0, hw = self.slant_inclined_path_coefficients(f, P + e)
+            h0, hw = self.slant_inclined_path_equivalent_height(f, P + e)
 
             # Use the zenit water-vapour method if the values of V_t
             # and h are provided
@@ -769,7 +774,7 @@ def gaseous_attenuation_inclined_path(
             gammaw = 0
 
         e = rho * T / 216.7
-        h0, hw = self.slant_inclined_path_coefficients(f, P + e)
+        h0, hw = self.slant_inclined_path_equivalent_height(f, P + e)
 
         if 5 < el and el < 90:
             h0_p = h0 * (np.exp(-h1 / h0) - np.exp(-h2 / h0))
@@ -861,8 +866,9 @@ def zenit_water_vapour_attenuation(self, *args, **kwargs):
     def gaseous_attenuation_approximation(self, *args, **kwargs):
         return _ITU676_10.gaseous_attenuation_approximation(*args, **kwargs)
 
-    def slant_inclined_path_coefficients(self, *args, **kwargs):
-        return _ITU676_10.slant_inclined_path_coefficients(*args, **kwargs)
+    def slant_inclined_path_equivalent_height(self, *args, **kwargs):
+        return _ITU676_10.slant_inclined_path_equivalent_height(*args,
+                                                                **kwargs)
 
     def gaseous_attenuation_terrestrial_path(self, *args, **kwargs):
         return _ITU676_10.gaseous_attenuation_terrestrial_path(*args, **kwargs)
@@ -1091,6 +1097,35 @@ def gaseous_attenuation_inclined_path(f, el, rho, P, T, h1, h2, mode='approx'):
     return prepare_output_array(val, type_output) * u.dB
 
 
+def slant_inclined_path_equivalent_height(f, p):
+    """ Computes the equivalent height to be used for oxygen and water vapour
+    gaseous attenuation computations.
+
+    Parameters
+    ----------
+    f : number or Quantity
+        Frequency (GHz)
+    p : number
+        Percentage of the time the gaseous attenuation value is exceeded.
+
+    Returns
+    -------
+    ho, hw : Quantity
+        Equivalent height for oxygen and water vapour (m)
+
+    References
+    --------
+    [1] Attenuation by atmospheric gases:
+    https://www.itu.int/rec/R-REC-P.676/en
+
+    """
+    type_output = type(f)
+    f = prepare_quantity(f, u.GHz, 'Frequency')
+
+    val = __model.slant_inclined_path_equivalent_height(f, p)
+    return prepare_output_array(val, type_output) * u.m
+
+
 @memory.cache
 def zenit_water_vapour_attenuation(lat, lon, p, f, V_t=None, h=None):
     """
diff --git a/itur/utils.py b/itur/utils.py
index faee0be..fd74b6f 100644
--- a/itur/utils.py
+++ b/itur/utils.py
@@ -329,9 +329,43 @@ def elevation_angle(h, lat_s, lon_s, lat_grid, lon_grid):
 
 
 def plot_in_map(data, lat=None, lon=None, lat_min=None, lat_max=None,
-                lon_min=None, lon_max=None, cbar_text='',
+                lon_min=None, lon_max=None, cbar_text='', ax=None,
                 **kwargs):
+    '''
+    Displays the value sin data in a map.
+
+    Either {lat, lon} or {lat_min, lat_max, lon_min, lon_max} need to be
+    provided as inputs.
+
+    Parameters
+    ----------
+    data : np.ndarray
+        Data values to be plotted.
+    lat : np.ndarray
+        Matrix with the latitudes for each point in data (deg N)
+    lon : np.ndarray
+        Matrix with the longitudes for each point in data (deg E)
+    lat_min :  float
+        Minimum latitude of the data (deg N)
+    lat_max :  float
+        Maximum latitude of the data (deg N)
+    lon_min :  float
+        Minimum longitude of the data (deg E)
+    lat_max :  float
+        Maximum longitude of the data (deg E)
+    cbar_text : string
+        Colorbar text caption.
+    ax : Axes
+        matplotlib axes where the data will be plotted.
+    **kwargs: dict
+        Key-value arguments that will be passed to the imshow function.
 
+
+    Returns
+    -------
+    m : Basemap
+        The map object generated by Basemap
+    '''
     import matplotlib.pyplot as plt
 
     try:
@@ -353,13 +387,15 @@ def plot_in_map(data, lat=None, lon=None, lat_min=None, lat_max=None,
         lon_max = np.max(lon)
         lon_min = np.min(lon)
 
-    ax = plt.subplot(111)
+    if ax is None:
+        ax = plt.subplot(111)
+
     m = Basemap(ax=ax, projection='cyl', llcrnrlat=lat_min,
                 urcrnrlat=lat_max, llcrnrlon=lon_min, urcrnrlon=lon_max,
                 resolution='l')
 
-    m.drawcoastlines(color='white', linewidth=0.5)
-    m.drawcountries(color='grey', linewidth=0.3)
+    m.drawcoastlines(color='grey', linewidth=0.8)
+    m.drawcountries(color='grey', linewidth=0.8)
     parallels = np.arange(-80, 81, 20)
     m.drawparallels(parallels, labels=[1, 0, 0, 1], dashes=[2, 1],
                     linewidth=0.2, color='white')
diff --git a/setup.py b/setup.py
index bff2a95..6aae49e 100644
--- a/setup.py
+++ b/setup.py
@@ -40,7 +40,7 @@
     # For a discussion on single-sourcing the version across setup.py and the
     # project code, see
     # https://packaging.python.org/en/latest/single_source_version.html
-    version='0.1.6',
+    version='0.2.0',
 
     # This is a one-line description or tagline of what your project does. This
     # corresponds to the "Summary" metadata field: