solve a major bug in the active mode on the mode coefficients. It was…

… introduced after the GMD paper. This commit also integrate numerical improvement of the full geometrical optics code.

smrt/core/lib.py

@@ -65,53 +65,64 @@ def len_atleast_1d(x):
 class smrt_diag(object):
     """Scipy.sparse is very slow for diagonal matrix and numpy has no good support for linear algebra. This diag class
     implements simple diagonal object without numpy subclassing (but without much features).
-    implements simple diagional object without numpy subclassing and without much features.
     It seems that proper subclassing numpy and overloading matmul is a very difficult problem."""
 
     __array_ufunc__ = None
 
     def __init__(self, arr):
         self.diag = np.atleast_1d(arr)
 
+        # self.shape = shape if shape is not None else (len(self.diag), len(self.diag))
+
     def diagonal(self):
         return self.diag
 
-    # @property
-    # def shape(self):
-    #     return self.values.shape
+    @property
+    def shape(self):
+        return (len(self.diag), len(self.diag))
 
     def __len__(self):
         return len(self.diag)
 
     def __rmatmul__(self, other):
         self.check_type(other)
+        # assert other.shape[1] == self.shape[0]
         return other * self.diag[np.newaxis, :]
 
     def __matmul__(self, other):
+        # assert self.shape[1] == other.shape[0]
         if isinstance(other, smrt_diag):
             # return a diagonal object
-            return smrt_diag(other.diag * self.diag)
+            return smrt_diag(other.diag * self.diag)  # , shape=(self.shape[0], other.shape[1]))
         else:
             self.check_type(other)
             # other must be an ndarray
             return other * self.diag[:, np.newaxis]
 
     def __rmul__(self, other):
         assert np.isscalar(other)
-        return other * self.diag
+        return smrt_diag(other * self.diag)
 
     def __mul__(self, other):
         assert np.isscalar(other)
-        return self.diag * other
+        return smrt_diag(self.diag * other)
 
     def __add__(self, other):
         if isinstance(other, smrt_diag):
             return smrt_diag(other.diag + self.diag)
         elif np.isscalar(other) and other == 0:
             return self
+        elif isinstance(other, np.ndarray):
+            assert other.shape == self.shape  # we do not allow broadcasting (not yet)...
+            other = other.copy()
+            other[np.diag_indices_from[other]] += self.diag
+            return other
         else:
             raise NotImplementedError
 
+    def __radd__(self, other):
+        return self.__add__(other)
+
     def __sub__(self, other):
         if isinstance(other, smrt_diag):
             return smrt_diag(other.diag - self.diag)
@@ -262,7 +273,6 @@ def compress(self, mode=None, auto_reduce_npol=False):
         else:
             raise NotImplementedError
 
-
     def __rmul__(self, other):
         return smrt_matrix(other * self.values)
 
@@ -352,7 +362,7 @@ def abs2(c):
     return c.real**2 + c.imag**2
 
 
-def generic_ft_even_matrix(phase_function, m_max):
+def generic_ft_even_matrix(phase_function, m_max, nsamples=None):
     """ Calculation of the Fourier decomposed of the phase or reflection or transmission matrix provided by the function.
 
     This method calculates the Fourier decomposition modes and return the output.
@@ -376,7 +386,10 @@ def generic_ft_even_matrix(phase_function, m_max):
     """
 
     # samples of dphi for fourier decomposition. Highest efficiency for 2^n. 2^2 ok
-    nsamples = 2**np.ceil(3 + np.log(m_max + 1) / np.log(2))
+    if nsamples is None:
+        nsamples = 2**np.ceil(3 + np.log(m_max + 1) / np.log(2))
+
+    assert nsamples > 2 * m_max
 
     # dphi must be evenly spaced from 0 to 2 * np.pi (but not including period), but we can use the symmetry of the phase function
     # to reduce the computation to 0 to pi (including 0 and pi) and mirroring for pi to 2*pi (excluding both)
@@ -389,8 +402,9 @@ def generic_ft_even_matrix(phase_function, m_max):
     npol = p.npol
 
     # mirror the phase function
+    assert len(p.values.shape) == 5
     p_mirror = p.values[:, :, -2:0:-1, :, :].copy()
-    if npol >= 3 :
+    if npol >= 3:
         p_mirror[0:2, 2] = -p_mirror[0:2, 2]
         p_mirror[2, 0:2] = -p_mirror[2, 0:2]
 
@@ -408,6 +422,7 @@ def generic_ft_even_matrix(phase_function, m_max):
 
     # m>=1 modes
     if npol == 2:
+        # the factor 2 comes from the change exp -> cos, i.e. exp(-ix) + exp(+ix)= 2 cos(x)
         ft_even_p[:, :, 1:] = ft_p[:, :, 1:m_max + 1].real * (2.0 / nsamples)
 
     else:

smrt/core/sensor.py

@@ -138,7 +138,9 @@ def active(frequency, theta_inc, theta=None, phi=None, polarization_inc=None, po
     if polarization_inc is None:
         polarization_inc = ['V', 'H']
 
-    sensor = Sensor(frequency, theta_inc, theta, phi, polarization_inc, polarization, channel_map=channel_map, name=name)
+    sensor = Sensor(frequency, theta_inc_deg=theta_inc, theta_deg=theta, phi_deg=phi,
+                    polarization_inc=polarization_inc, polarization=polarization,
+                    channel_map=channel_map, name=name)
 
     sensor.basic_checks()
 

smrt/interface/coherent_flat.py

@@ -10,6 +10,7 @@
 from smrt.core.globalconstants import C_SPEED
 from smrt.core.lib import smrt_matrix, abs2
 from smrt.core.fresnel import fresnel_coefficients
+from smrt.core.error import SMRTError
 
 
 def process_coherent_layers(snowpack, emmodel_list, sensor):
@@ -21,7 +22,7 @@ def process_coherent_layers(snowpack, emmodel_list, sensor):
 
     if not np.any(coherent_layers):
         return snowpack, emmodel_list
-    
+
     snowpack = snowpack.copy()
 
     if coherent_layers[-1]:
@@ -31,16 +32,16 @@ def process_coherent_layers(snowpack, emmodel_list, sensor):
 
     for l in np.flatnonzero(coherent_layers[:-1])[::-1]:  # reverse the processing to safely delete the snowpack layer and interface
         print("coherent layer:", l)
-        if coherent_layers[l-1]:
+        if coherent_layers[l - 1]:
             raise SMRTError("Two sucessive layers are coherent, this is not yet supported")
         # create a coherent interface
-        coherent_interface = CoherentFlat(snowpack.interfaces[l:l+2], snowpack.layers[l], effective_permittivity[l])
+        coherent_interface = CoherentFlat(snowpack.interfaces[l:l + 2], snowpack.layers[l], effective_permittivity[l])
         # set the next interface to coherent
         snowpack.interfaces[l + 1] = coherent_interface
         # delete the layer to be deleted
         snowpack.delete(l)  # delete layer and interface l
         emmodel_list.pop(l)
-           
+
     return snowpack, emmodel_list
 
 
@@ -78,7 +79,7 @@ def specular_reflection_matrix(self, frequency, eps_1, eps_2, mu1, npol):
         reflection_coefficients[1] = abs2(R_h)
 
         if npol >= 3:
-            reflection_coefficients[2] = (R_v*np.conj(R_h)).real   # TsangI  Eq 7.2.93
+            reflection_coefficients[2] = (R_v * np.conj(R_h)).real   # TsangI  Eq 7.2.93
 
         return reflection_coefficients
 
@@ -105,23 +106,22 @@ def coherent_transmission_matrix(self, frequency, eps_1, eps_2, mu1, npol):
 
         transmission_coefficients = smrt_matrix.ones((npol, len(mu1)))
 
-        nt = np.sqrt(eps_2/eps_1).real
-        transmission_coefficients[0] = abs2(T_v) * mu_t/mu1 / nt  # for the coef see TsangIII 2.1.140b
-        transmission_coefficients[1] = abs2(T_h) * mu_t/mu1 * nt  # for the coef see TsangIII 2.1.140a
+        nt = np.sqrt(eps_2 / eps_1).real
+        transmission_coefficients[0] = abs2(T_v) * mu_t / mu1 / nt  # for the coef see TsangIII 2.1.140b
+        transmission_coefficients[1] = abs2(T_h) * mu_t / mu1 * nt  # for the coef see TsangIII 2.1.140a
 
         if npol >= 3:
             # this part is to be confirmed.
             R_v = (R01_v + R1t_v * exp_2kd) / (1 + R01_v * R1t_v * exp_2kd)
             R_h = (R01_h + R1t_h * exp_2kd) / (1 + R01_h * R1t_h * exp_2kd)
 
-            transmission_coefficients[2] = mu_t / mu1 * ((1+R_v)*np.conj(1+R_h)).real  # TsangI  Eq 7.2.95
+            transmission_coefficients[2] = mu_t / mu1 * ((1 + R_v) * np.conj(1 + R_h)).real  # TsangI  Eq 7.2.95
 
         return transmission_coefficients
 
-
     def _prepare_computation(self, frequency, eps_1, eps_2, mu1):
 
-            # convert to Tsang's notation. See TsangI, pages 207, Eq. 5.2.14
+        # convert to Tsang's notation. See TsangI, pages 207, Eq. 5.2.14
         eps_0 = eps_1
         eps_1 = self.permittivity
         eps_t = eps_2
@@ -133,13 +133,16 @@ def _prepare_computation(self, frequency, eps_1, eps_2, mu1):
         k_1 = 2 * np.pi / C_SPEED * frequency * np.sqrt(eps_1)
 
         phase = k_1 * mu_1 * self.layer.thickness
-        assert np.all(phase.imag >=0)
+        assert np.all(phase.imag >= 0)
 
         incoherent = phase.real > 3 * np.pi / 4  # we consider coherency up to 3 pi / 2 like in MEMLS
 
         phase[incoherent].real = 0
 
-        exp_kd = np.exp(1j  * phase)
-        exp_2kd = np.exp(2j  * phase)
+        exp_kd = np.exp(1j * phase)
+        exp_2kd = np.exp(2j * phase)
 
-        return R01_v, R01_h, R1t_v, R1t_h, exp_kd, exp_2kd, mu_t
+        return R01_v, R01_h, R1t_v, R1t_h, exp_kd, exp_2kd, mu_t
+
+    def diffuse_transmission_matrix(self, frequency, eps_1, eps_2, mu_s, mu_i, dphi, npol):
+        return smrt_matrix(0)

smrt/interface/flat.py

@@ -47,3 +47,5 @@ def coherent_transmission_matrix(self, frequency, eps_1, eps_2, mu1, npol):
 
         return fresnel_transmission_matrix(eps_1, eps_2, mu1, npol)
 
+    def diffuse_transmission_matrix(self, frequency, eps_1, eps_2, mu_s, mu_i, dphi, npol):
+        return smrt_matrix(0)