diff --git a/smrt/interface/iiem_fung2002.py b/smrt/interface/iiem_fung2002.py
index 914634a..924cd67 100644
--- a/smrt/interface/iiem_fung2002.py
+++ b/smrt/interface/iiem_fung2002.py
@@ -219,16 +219,18 @@ def diffuse_reflection_matrix(self, frequency, eps_1, eps_2, mu_s, mu_i, dphi, n
             Rv, Rh, _ = fresnel_coefficients(eps_1, eps_2, mu_i)
             Rvh = (Rv - Rh) / 2
             ks2 = ks**2
-            svh = self.double_integral(k, ks2, mu_i, eps_2, Rvh, n, self.N_integral)
+
+            # only compute for incident angle
+            # reshape to match angles shape set initially
+            svh = self.double_integral(k, ks2, mu_i, eps_2, Rvh, n, self.N_integral).reshape(1, 1, mu_i.shape[2])
             # reshape to match angles shape set initially
-            svh = svh.reshape(1, 1, mu_i.shape[2], 1)
+            mu_i = mu_i.squeeze(axis = -1)
+            sin_i = sin_i.squeeze(axis = -1)
 
             if self.shadow_correction:
                 s = 1 / (1 + shadow_function(mean_square_slope, mu_i / sin_i) * 2)
                 svh *= s
 
-            # reshape again to match result final shape
-            svh = svh.reshape(1, 1, mu_i.shape[2])
             reflection_coefficients[0, 1] = svh / (4 * np.pi * mu_i)
             reflection_coefficients[1, 0] = svh / (4 * np.pi * mu_i)
 
diff --git a/smrt/rtsolver/iterative_second_order.py b/smrt/rtsolver/iterative_second_order.py
index f684ac0..c18a0ce 100644
--- a/smrt/rtsolver/iterative_second_order.py
+++ b/smrt/rtsolver/iterative_second_order.py
@@ -496,12 +496,12 @@ def compute_scattering_layer_ground(
         n_mu_i = len(mus_i)
 
         phase_mu_int_mu = emmodel.ft_even_phase(mu_int_sym, mu_i_sym, m_max) / (4 * np.pi)
-        # phase_mu_mu_int = emmodel.ft_even_phase(mu_i_sym, mu_int_sym, m_max) / (4 * np.pi)
+        phase_mu_mu_int = emmodel.ft_even_phase(mu_i_sym, mu_int_sym, m_max) / (4 * np.pi)
 
         R1 = Rbottom_diff_int["i_int"][:, :, :, n_mu_i:, n_stream:]  # R(mu_i, mu_int)
         P1 = phase_mu_int_mu[:, :, :, 0:n_stream, 0:n_mu_i]  # P(-mu_int, -mu_i)
-        R2 = Rbottom_diff_int["i_int"][:, :, :, n_mu_i:, 0:n_stream]  # R(mu_i, -mu_int)
-        P2 = phase_mu_int_mu[:, :, :, n_stream:, 0:n_mu_i]  # P(mu_int, -mu_i)
+        R2 = Rbottom_diff_int["int_i"][:, :, :, n_stream:, 0:n_mu_i]  # R(mu_int, -mu_i)
+        P2 = phase_mu_mu_int[:, :, :, n_mu_i:, n_stream:]  # P(mu_i, mu_int)
 
         sum_e = 0
         for mu, w, i in zip(mu_int, weight, range(n_stream)):
@@ -513,7 +513,7 @@ def compute_scattering_layer_ground(
             E = compute_E(mus_i, mu, ke, layer_optical_depth, layer_optical_depth_ln_ground)
             sum_e += w * (
                 (E * compute_integral_phi(R1[:, :, :, :, i], P1[:, :, :, i, :], m_max, len_mu, npol, np.pi))
-                + (E * compute_integral_phi(R2[:, :, :, :, i], P2[:, :, :, i, :], m_max, len_mu, npol, np.pi))
+                + (E * compute_integral_phi(R2[:, :, :, i, :], P2[:, :, :, :, i], m_max, len_mu, npol, np.pi))
             )
 
         return sum_e @ I_l