update the tests after the change of Fresnel coefficients

CHANGELOG.md

@@ -8,4 +8,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 ### Added
 	- volumetric_liquid_water is a new (and recommended) way to set the amount of water in a layer using make_snowpack.
 
-	- the 'emmodel' argument in make_model can now be a dict mapping a different emmodel for each sort of layer medium. Useful for snow + sea-ice for instance when the emmodels must be different for snow and ice.
+	- the 'emmodel' argument in make_model can now be a dict mapping of different emmodels for each sort of layer medium. Useful for snow + sea-ice for instance when the emmodels must be different for snow and ice.
+
+	- a new function in make_medium.py to create a water body (lake or open ocean): from smrt import make_water_body  
+
+
+### Changed
+	- Fresnel coefficients for the reflection and transmission on flat interface are now calculated with a more rigorous equation for (very) lossly materials. This could affect some simulations but with very little effect 60°.
+

smrt/rtsolver/dort.py

@@ -699,6 +699,8 @@ def solve(self, m, compute_coherent_only, debug_A=False):
                     print(np.any(mask, axis=1))
                     print(beta[np.any(mask, axis=1)])
                     print(beta)
+                else:
+                    E = E.real
 
             beta = beta.real
 

smrt/test/test_coherent_layer.py

@@ -37,7 +37,7 @@ def test_snowpack_with_coherent_layer():
     #assert abs(res.TbH() - 196.83495992559307) < 1e-4
 
     # the new values come form the correction of 917->916.7
-    assert abs(res.TbV() - 261.06421847662216) < 1e-4
-    assert abs(res.TbH() - 196.8660443556061) < 1e-4
+    assert abs(res.TbV() - 261.0633483757312) < 1e-4
+    assert abs(res.TbH() - 196.8659636937278) < 1e-4
 
 

smrt/test/test_iba_sea_ice.py

@@ -7,8 +7,8 @@
 from smrt.inputs.make_medium import make_ice_column, bulk_ice_density
 from smrt.core.interface import Interface
 
-#test if this configuration gives values as originally produced by examples/iba_sea_ice.py
-#same structure as test_integration_iba.py
+# test if this configuration gives values as originally produced by examples/iba_sea_ice.py
+# same structure as test_integration_iba.py
 
 
 def setup_seaice():
@@ -17,7 +17,8 @@ def setup_seaice():
     n_max_stream = 64
 
     thickness = np.array([1.5 / l] * l)  # ice is 1.5m thick
-    temperature = np.linspace(273.15 - 20., 273.15 - 1.8, l)  # temperature gradient in the ice from -20 deg C at top to freezing temperature of water at bottom (-1.8 deg C)
+    # temperature gradient in the ice from -20 deg C at top to freezing temperature of water at bottom (-1.8 deg C)
+    temperature = np.linspace(273.15 - 20., 273.15 - 1.8, l)
     salinity = np.linspace(2., 10., l) * PSU  # salinity profile ranging from salinity=2 at the top to salinity=10 at the bottom of the ice
 
     return l, n_max_stream, thickness, temperature, salinity
@@ -34,26 +35,25 @@ def test_oneconfig_for_firstyear_sea_ice():
                                  thickness=thickness,
                                  temperature=temperature,
                                  microstructure_model="exponential",
-                                 brine_inclusion_shape="spheres", #inclusion_shape can be "spheres" or "random_needles", or "mix_spheres_needles"
-                                 salinity=salinity, #either 'salinity' or 'brine_volume_fraction' should be given for sea ice; if salinity is given, brine volume fraction is calculated in the model; if none is given, ice is treated as fresh water ice
+                                 brine_inclusion_shape="spheres",  # inclusion_shape can be "spheres" or "random_needles", or "mix_spheres_needles"
+                                 salinity=salinity,  # either 'salinity' or 'brine_volume_fraction' should be given for sea ice; if salinity is given, brine volume fraction is calculated in the model; if none is given, ice is treated as fresh water ice
                                  corr_length=p_ex,
                                  add_water_substrate="ocean"
                                  )
 
     # create the sensor
     sensor = sensor_list.passive(1.4e9, 40.)
-    n_max_stream= 128
-    m = make_model("iba", "dort", rtsolver_options ={"n_max_stream": n_max_stream})
+    n_max_stream = 128
+    m = make_model("iba", "dort", rtsolver_options={"n_max_stream": n_max_stream})
 
     # run the model
     res = m.run(sensor, ice_column)
 
     print(res.TbV(), res.TbH())
-    #absorption with effective permittivity
-    assert abs(res.TbV() - 256.01165061598317) < 1e-4
-    assert abs(res.TbH() - 228.47447378338745) < 1e-4
+    # absorption with effective permittivity
+    assert abs(res.TbV() - 256.0170296269674) < 1e-4
+    assert abs(res.TbH() - 228.4566040823167) < 1e-4
 
-
 
 def test_oneconfig_for_multiyear_sea_ice():
     # prepare inputs
@@ -76,16 +76,16 @@ def test_oneconfig_for_multiyear_sea_ice():
     # create the sensor
     sensor = sensor_list.passive(1.4e9, 40.)
 
-    n_max_stream= 128
-    m = make_model("iba", "dort", rtsolver_options ={"n_max_stream": n_max_stream})
+    n_max_stream = 128
+    m = make_model("iba", "dort", rtsolver_options={"n_max_stream": n_max_stream})
 
     # run the model
     res = m.run(sensor, ice_column)
 
     print(res.TbV(), res.TbH())
     # absorption with effective permittivity
-    assert abs(res.TbV() - 257.5689162188296) < 1e-4
-    assert abs(res.TbH() - 232.03924683304172) < 1e-4
+    assert abs(res.TbV() - 257.57209000420636) < 1e-4
+    assert abs(res.TbH() - 232.01555447145563) < 1e-4
 
 
 def test_equivalence_porosity_density():
@@ -111,7 +111,7 @@ def test_equivalence_porosity_density():
                                   )
 
     # Same, but giving the density instead:
-    density = [bulk_ice_density(temp, salt , porosity) for temp, salt in zip(temperature, salinity)]
+    density = [bulk_ice_density(temp, salt, porosity) for temp, salt in zip(temperature, salinity)]
 
     ice_column2 = make_ice_column(ice_type,
                                   thickness=thickness,

smrt/test/test_integration_iba.py

@@ -9,8 +9,8 @@
 # Ghi: rapid hack, should be splitted in different functions
 #
 
-def setup_snowpack(l):
 
+def setup_snowpack(l):
 
     nl = l//2  # // Forces integer division
     thickness = np.array([0.1, 0.1]*nl)
@@ -25,7 +25,8 @@ def setup_snowpack(l):
                              density=density,
                              temperature=temperature,
                              corr_length=p_ex)
-    return  snowpack
+    return snowpack
+
 
 def test_iba_oneconfig_passive():
 
@@ -42,13 +43,12 @@ def test_iba_oneconfig_passive():
     res = m.run(radiometer, snowpack)
 
     print(res.TbV(), res.TbH())
-    #absorption with effective permittivity
+    # absorption with effective permittivity
     # abs(res.TbV() - 248.08794944809972) < 1e-4
     # abs(res.TbH() - 237.3056263719142) < 1e-4
 
-    assert abs(res.TbV() - 248.08744066791073) < 1e-4
-    assert abs(res.TbH() - 237.30720491883298) < 1e-4
-
+    assert abs(res.TbV() - 248.08374547409588) < 1e-4
+    assert abs(res.TbH() - 237.30435496083572) < 1e-4
 
 
 def test_iba_oneconfig_active():
@@ -67,7 +67,7 @@ def test_iba_oneconfig_active():
 
     print(res.sigmaVV_dB(), res.sigmaHH_dB(), res.sigmaHV_dB())
 
-    assert abs(res.sigmaVV_dB() - (-24.04497237)) < 1e-3
-    assert abs(res.sigmaHH_dB() - (-24.41628343)) < 1e-3
-    assert abs(res.sigmaHV_dB() - (-51.53673914)) < 1e-3
+    assert abs(res.sigmaVV_dB() - (-24.044882546524693)) < 1e-3
+    assert abs(res.sigmaHH_dB() - (-24.416295329469907)) < 1e-3
+    assert abs(res.sigmaHV_dB() - ( -51.544272924876886)) < 1e-3
 

smrt/test/test_integration_iba_original.py

@@ -42,5 +42,5 @@ def test_iba_oneconfig():
     print(res.TbV(), res.TbH())
 
     #original absorption (Maetzler 1998)
-    assert abs(res.TbV() - 247.92344524715216) < 1e-4
-    assert abs(res.TbH() - 237.0863426896562) < 1e-4
+    assert abs(res.TbV() - 247.91935638633657) < 1e-4
+    assert abs(res.TbH() - 237.08319097431894) < 1e-4

smrt/test/test_physics_law.py

@@ -35,12 +35,13 @@ def _do_test_isothermal_universe(pc, thickness_1):
 
     T = 265
 
-    substrate = make_soil('soil_wegmuller', permittivity_model=complex(10,1), roughness_rms=0.001, temperature=T)
+    substrate = make_soil('soil_wegmuller', permittivity_model=complex(10, 1), roughness_rms=0.001, temperature=T)
 
     atmosphere = SimpleIsotropicAtmosphere(tbdown=T, tbup=0, trans=1)
 
     snowpack = make_snowpack([0.3, thickness_1], "exponential",
                              density=[200, 300], temperature=T, corr_length=pc,
+                             ice_permittivity_model=complex(1.7, 0.00001),
                              substrate=substrate, atmosphere=atmosphere)
 
     # create the sensor
@@ -59,14 +60,15 @@ def _do_test_isothermal_universe(pc, thickness_1):
 
 def _do_test_kirchoff_law(pc, thickness_1):
 
-    T = 265
+    T = 265.
 
     substrate = make_soil('soil_wegmuller', permittivity_model=complex(10, 1), roughness_rms=0.001, temperature=T)
 
     atmosphere1K = SimpleIsotropicAtmosphere(tbdown=1, tbup=0, trans=1)
 
     snowpack = make_snowpack([0.3, thickness_1], "exponential",
                              density=[200, 300], temperature=T, corr_length=pc,
+                             ice_permittivity_model=complex(1.7, 0.00001),
                              substrate=substrate)
 
     # create the sensor

smrt/test/test_soil.py

@@ -54,8 +54,8 @@ def test_soil_wegmuller_dobson85():
     res = run_model(snowpack)
 
     print(res.TbV(), res.TbH())
-    assert abs(res.TbV() - 262.60009290172155) < 1e-4
-    assert abs(res.TbH() - 255.8655605977706) < 1e-4
+    assert abs(res.TbV() - 262.55457107119486) < 1e-4
+    assert abs(res.TbH() - 255.81725907587176) < 1e-4
     # note value from DMRTML Fortran running in the same conditions:
     # H=255.88187817295605 V=262.60345275739024
 
@@ -72,7 +72,7 @@ def test_soil_wegmuller_montpetit2008():
     res = run_model(snowpack)
 
     print(res.TbV(), res.TbH())
-    assert abs(res.TbV() - 262.45644269568135) < 1e-4
-    assert abs(res.TbH() - 255.7131655561391) < 1e-4
+    assert abs(res.TbV() - 262.4543081568107) < 1e-4
+    assert abs(res.TbH() - 255.71089039573724) < 1e-4
     # note value from DMRTML Fortran running in the same conditions:
     # H=255.88187817295605 V=262.60345275739024