<jupyter_start><jupyter_text>Sensitivity analysis============__Goal__: - run sensitivity analysis to show the impact of a given parameter on the SMRT output  __Learning__:  SMRT is able to iterate on several arguments when it is unambiguous. For instance, a sensor with multiple frequencies, angles or polarizations is automatically understood. The `result` contains all the values which can be accessed with arguments in TbV() and similar functions. E.g. TbV(frequency=37e9)This is similar when a list of snowpacks is given to `run`. The `result` contains all the computations. The 'snowpack' dimension is automatically added but we can also propose a custom name for this dimension.In the following, we show different approaches to conduct sensitivity studies that  you can run and then apply to a study case of your choice: - take the Dome C snowpack and study the sensitivity of TbH 55° to superficial density - take any snowpack previously defined and investigated the sensivitiy to liquid_water - etc<jupyter_code>import numpy as np

import matplotlib.pyplot as plt
%matplotlib notebook

from smrt import make_model, make_snowpack, sensor_list<jupyter_output><empty_output><jupyter_text>Build a list of snowpack--------------------------------The key idea is to build a list of snowpack. E.g. we want to test the sensitivity of TB's to the radius. We first build a list of snowpack with different radius.<jupyter_code># prepare the snowpack
density = 300.0
radius = np.arange(0.05, 0.5, 0.01) *1e-3  # from 0 to 0.5mm

# the NAIVE APPROACH:

snowpack = list()
for x in radius:
    sp = make_snowpack([1000.0], "sticky_hard_spheres", 
                       density=density, temperature=265, 
                       radius=x, stickiness=0.15)
    snowpack.append(sp)<jupyter_output><empty_output><jupyter_text>In simple cases (as this one), it is easier to use "list comprehension", a nice python feature to create list.<jupyter_code># with list comprehension
snowpack = [make_snowpack([1000.0], "sticky_hard_spheres", 
                          density=density, temperature=265,
                          radius=x, stickiness=0.15) for x in radius]
# prepare the sensor and model

model = make_model("iba", "dort")
sensor = sensor_list.passive(37e9, 55)

#run!<jupyter_output><empty_output><jupyter_text>Now we have a list of snowpacks, we want to call the model for each snowpack. We can use list comprehension again.<jupyter_code>results = [model.run(sensor, sp) for sp in snowpack]<jupyter_output><empty_output><jupyter_text>This return a list of results. To extract the TB V for each result can be done with another list comprehension. And then we plot the results.<jupyter_code>tbv = [res.TbV() for res in results]
plt.figure()
plt.plot(radius, tbv)<jupyter_output><empty_output><jupyter_text>Nice ? We can do much better because `Model` can directly run on a list of snowpacks. It does not return a list of results, but **a unique result with a new coordinate** which is much more convenient.<jupyter_code>results = model.run(sensor, snowpack, snowpack_dimension=('radius', radius))
print(type(results))  # look results is a Result, not a list
print(results.coords) # look, we have several coordinates, one is call corr_legn<jupyter_output><empty_output><jupyter_text>This is more compact and nicer, `results` explicitly show the radius dimension. Plotting is thus easier:<jupyter_code>plt.figure()
plt.plot(results.radius, results.TbV())<jupyter_output><empty_output><jupyter_text>And it is easy to save the result to disk:<jupyter_code>results.save("radius-sensitivity.nc")<jupyter_output><empty_output><jupyter_text>Recap:---------<jupyter_code>snowpack = [make_snowpack([1000.0], "sticky_hard_spheres", density=density, temperature=265, radius=x, stickiness=0.15) for x in radius]

model = make_model("iba", "dort")
sensor = sensor_list.passive([19e9, 37e9], 55)

results = model.run(sensor, snowpack, snowpack_dimension=('radius', radius))

plt.figure()
plt.plot(results.radius, results.TbV(frequency=19e9), label="19 GHz")
plt.plot(results.radius, results.TbV(frequency=37e9), label="37 GHz")
plt.legend()
results.TbV()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup