implementing geometrical optics for rough interfaces within ice column about smrt HOT 4 OPEN

# import libraries
import numpy as np

# import smrt and related functions
from smrt import make_ice_column, make_snowpack, make_model, sensor
from smrt import PSU # g/kg -> kg/kg conversion
from smrt.atmosphere.simple_isotropic_atmosphere import SimpleIsotropicAtmosphere
from smrt.core.interface import make_interface # import lib for roughness

# ice inputs
l           = 2                             # number ice layers
thickness   = np.array([.5,.5])             # ice thickness in m
radius      = np.array([.0001,.0001])       # particle radius
stickiness  = np.array([0.3, 0.3])          # 'tau' 
temperature = np.array([260, 265])          # ice temperature in K
salinity    = np.array([5.1, 5.4])*PSU      # ice salinity profile in kg/kg
density     = np.array([0.924,0.924])       # ice density profile in kg/m^3
mean_square_slope = 0.1                       # MSS (unitless)

# atmosphere
atmos = SimpleIsotropicAtmosphere(10., 3., 0.90)

# create a first-year sea ice column:
ice_type = 'firstyear' # first-year or multi-year sea ice

rough_interface = make_interface("geometrical_optics",\
                                 mean_square_slope=mean_square_slope)

ice_column = make_ice_column(ice_type              = ice_type,
                             thickness             = thickness, # meters
                             temperature           = temperature,
                             microstructure_model  = "sticky_hard_spheres",
                             brine_inclusion_shape = "spheres", 
                             salinity              = salinity, 
                             brine_volume_fraction = .02,
                             radius                = radius,
                             stickiness            = stickiness,
                             density               = density,
                             add_water_substrate   = "ocean" #see comment below
                            )
ice_column.interfaces[0]=rough_interface
#ice_column.interfaces[1]=rough_interface

# snow inputs:
l_s           = 2                            # number of layers
thickness_s   = np.array([.12,.12])          # thickness per layer
density_s     = np.array([216,100])          # density profile in kg/m^3
radius_s      = np.array([.0004, .00025])    # particle radius
stickiness_s  = np.array([0.20, 0.20])       # 'tau' (stickiness)
temperature_s = np.array([253, 252])         # temperature

# create the snowpack
snowpack = make_snowpack(thickness            = thickness_s,
                         microstructure_model = "sticky_hard_spheres",
                         density              = density_s,
                         temperature          = temperature_s,
                         radius               = radius_s,
                         stickiness           = stickiness_s)

snowpack.atmosphere = atmos# from test

#add snowpack on top of ice column:
medium = snowpack + ice_column

# create geometric and EM parameters
theta        = 55           # Earth incidence angle
freq         = 6E9     # range of frequencies  (Hz)
sensor=sensor.passive(freq, theta)

# make model
n_max_stream = 256           # TB calcis more accurate if # of streams increased
                            # (currently: default = 32); 
                            # needs to be increased when using > 1 snow layer 
m = make_model("iba", "dort", rtsolver_options ={"n_max_stream": n_max_stream},
               emmodel_options=dict(dense_snow_correction="auto"))

# run the model for bare sea ice:
res1 = m.run(sensor, ice_column)
# run the model for snow-covered sea ice:
res2 = m.run(sensor, medium)

# print TBs at horizontal and vertical polarization:
print(res1.TbH(), res1.TbV())
print(res2.TbH(), res2.TbV())

207.19735077300595 245.69465007105475
200.60817061769922 227.73977475803218

5427.083559428849 4236.1671434738155
-16228.632555060074 -12573.256466463112