2. Clearsky radiative transfer
Zeeman
1import pyarts
2import numpy as np
3import matplotlib.pyplot as plt
4
5ws = pyarts.workspace.Workspace()
6
7# %% Sampled frequency range
8
9line_f0 = 118750348044.712
10ws.frequency_grid = np.linspace(-50e6, 50e6, 1001) + line_f0
11
12# %% Species and line absorption
13
14ws.absorption_speciesSet(species=["O2-66"])
15ws.ReadCatalogData()
16ws.absorption_bandsSelectFrequency(fmin=40e9, fmax=120e9, by_line=1)
17ws.absorption_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
18ws.WignerInit()
19
20# %% Use the automatic agenda setter for propagation matrix calculations
21ws.propagation_matrix_agendaAuto()
22
23# %% Grids and planet
24
25ws.surface_fieldSetPlanetEllipsoid(option="Earth")
26ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0
27ws.atmospheric_fieldRead(
28 toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
29)
30ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37")
31
32# %% Checks and settings
33
34ws.spectral_radiance_unit = "Tb"
35ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground")
36ws.spectral_radiance_surface_agendaSet(option="Blackbody")
37
38# %% Core calculations
39
40pos = [100e3, 0, 0]
41los = [180.0, 0.0]
42ws.ray_pathGeometric(pos=pos, los=los, max_step=1000.0)
43ws.spectral_radianceClearskyEmission()
44ws.spectral_radianceApplyUnitFromSpectralRadiance()
45
46# %% Show results
47
48plt.plot((ws.frequency_grid - line_f0) / 1e6, ws.spectral_radiance + 0)
49plt.xlabel("Frequency offset [MHz]")
50plt.ylabel("Spectral radiance [K]")
51plt.title(f"Zeeman effect of {round(line_f0/1e6)} MHz O$_2$ line")
52
53# %% Test
54
55assert np.allclose(
56 ws.spectral_radiance[::100],
57 np.array(
58 [
59 [2.27784834e02, -2.23635319e-04, -3.77001634e-04, 5.69266632e-02],
60 [2.30863853e02, -3.48081919e-04, -5.82851478e-04, 7.04138002e-02],
61 [2.34806521e02, -6.21501805e-04, -1.02884266e-03, 9.33836080e-02],
62 [2.40360801e02, -1.41669913e-03, -2.28974401e-03, 1.40029968e-01],
63 [2.49782407e02, -5.60506178e-03, -8.33739163e-03, 2.69735122e-01],
64 [2.07248447e02, -4.30842579e00, -2.14475851e01, 1.25057382e-05],
65 [2.49781840e02, -5.60734272e-03, -8.34027514e-03, -2.69812575e-01],
66 [2.40359708e02, -1.41773814e-03, -2.29128954e-03, -1.40106546e-01],
67 [2.34804943e02, -6.22167771e-04, -1.02988483e-03, -9.34596090e-02],
68 [2.30861802e02, -3.48574780e-04, -5.83642167e-04, -7.04898114e-02],
69 [2.27782317e02, -2.24029336e-04, -3.77643130e-04, -5.70026253e-02],
70 ]
71 ),
72), "Values have drifted from expected results in spectral radiance"
Zeeman sensor
1import pyarts
2import numpy as np
3import matplotlib.pyplot as plt
4
5ws = pyarts.workspace.Workspace()
6
7# %% Sampled frequency range
8
9line_f0 = 118750348044.712
10ws.frequency_grid = np.linspace(-50e6, 50e6, 1001) + line_f0
11
12# %% Species and line absorption
13
14ws.absorption_speciesSet(species=["O2-66"])
15ws.ReadCatalogData()
16ws.absorption_bandsSelectFrequency(fmin=40e9, fmax=120e9, by_line=1)
17ws.absorption_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
18ws.WignerInit()
19
20# %% Use the automatic agenda setter for propagation matrix calculations
21ws.propagation_matrix_agendaAuto()
22
23# %% Grids and planet
24
25ws.surface_fieldSetPlanetEllipsoid(option="Earth")
26ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0
27ws.atmospheric_fieldRead(
28 toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
29)
30ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37")
31
32# %% Checks and settings
33
34ws.spectral_radiance_unit = "Tb"
35ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground")
36ws.spectral_radiance_surface_agendaSet(option="Blackbody")
37ws.ray_path_observer_agendaSet(option="Geometric")
38ws.spectral_radiance_observer_agendaSet(option="EmissionUnits")
39
40# %% Set up a sensor with Gaussian channel widths on individual frequency ranges
41
42pos = [100e3, 0, 0]
43los = [180.0, 0.0]
44ws.measurement_sensorSimpleGaussian(fwhm=1e5, pos=pos, los=los, pol="RC")
45
46# %% Core calculations
47
48result = pyarts.arts.Vector()
49result_jac = pyarts.arts.Matrix()
50ws.measurement_vectorFromSensor(result, result_jac)
51
52# %% Show results
53
54plt.plot((ws.frequency_grid - line_f0) / 1e6, result)
55plt.xlabel("Frequency offset [MHz]")
56plt.ylabel("Spectral radiance [K]")
57plt.title(
58 f"Zeeman effect of {round(line_f0/1e6)} MHz O$_2$ line with Gaussian channels on individual grids"
59)
60
61# %% Test
62
63assert np.allclose(
64 result[::100],
65 np.array(
66 [
67 227.78646795,
68 230.8638575,
69 234.80652899,
70 240.36081974,
71 249.78247057,
72 207.62113428,
73 249.78190355,
74 240.35972683,
75 234.80495168,
76 230.86180615,
77 227.78395156,
78 ]
79 ),
80)
Zeeman transmission
1import pyarts
2import numpy as np
3import matplotlib.pyplot as plt
4
5ws = pyarts.workspace.Workspace()
6
7# %% Sampled frequency range
8
9line_f0 = 118750348044.712
10nf = 1001
11ws.frequency_grid = np.linspace(-50e6, 50e6, nf) + line_f0
12
13# %% Species and line absorption
14
15ws.absorption_speciesSet(species=["O2-66"])
16ws.ReadCatalogData()
17ws.absorption_bandsSelectFrequency(fmin=40e9, fmax=120e9, by_line=1)
18ws.absorption_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
19ws.WignerInit()
20
21# %% Use the automatic agenda setter for propagation matrix calculations
22ws.propagation_matrix_agendaAuto()
23
24# %% Grids and planet
25
26ws.surface_fieldSetPlanetEllipsoid(option="Earth")
27ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0
28ws.atmospheric_fieldRead(
29 toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
30)
31ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37")
32
33# %% Checks and settings
34
35ws.spectral_radiance_unit = "1"
36ws.spectral_radiance_space_agendaSet(option="Transmission")
37ws.spectral_radiance_surface_agendaSet(option="Transmission")
38
39# %% Core calculations
40
41pos = [100e3, 0, 0]
42los = [180.0, 0.0]
43ws.ray_pathGeometric(pos=pos, los=los, max_step=1000.0)
44ws.spectral_radianceClearskyTransmission()
45
46# %% Show results
47
48plt.semilogy((ws.frequency_grid - line_f0) / 1e6, ws.spectral_radiance)
49plt.xlabel("Frequency offset [MHz]")
50plt.ylabel("Spectral radiance [K]")
51plt.title(f"Zeeman effect of {round(line_f0/1e6)} MHz O$_2$ line")
52
53# %% Test
54
55assert np.allclose(
56 ws.spectral_radiance[::100],
57 np.array(
58 [
59 [3.48074339e-06, 7.79626187e-11, 1.29462746e-10, -4.58316733e-08],
60 [1.71316764e-06, 5.95493071e-11, 9.78633101e-11, -2.77698065e-08],
61 [6.97243885e-07, 4.28547196e-11, 6.92129313e-11, -1.47794484e-08],
62 [2.01811291e-07, 2.79796924e-11, 4.36322508e-11, -6.24892024e-09],
63 [2.58887859e-08, 1.50722487e-11, 2.10895501e-11, -1.53914529e-09],
64 [3.25340224e-13, 4.13333797e-14, -3.14840689e-13, 8.41284772e-18],
65 [2.59107134e-08, 1.50971400e-11, 2.11231675e-11, 1.54394772e-09],
66 [2.02108653e-07, 2.80648334e-11, 4.37633416e-11, 6.28221755e-09],
67 [6.98656898e-07, 4.30420221e-11, 6.95130288e-11, 1.48859032e-08],
68 [1.71750402e-06, 5.98858575e-11, 9.84130185e-11, 2.80160699e-08],
69 [3.49118434e-06, 7.85004414e-11, 1.30351250e-10, 4.63074542e-08],
70 ]
71 ),
72), "Values have drifted from expected results in spectral radiance"
Zeeman sun
1import pyarts
2import numpy as np
3import matplotlib.pyplot as plt
4
5ws = pyarts.workspace.Workspace()
6
7# %% Sampled frequency range
8
9line_f0 = 118750348044.712
10ws.frequency_grid = [line_f0]
11
12# %% Species and line absorption
13
14ws.absorption_speciesSet(species=["O2-66"])
15ws.ReadCatalogData()
16ws.absorption_bandsSelectFrequency(fmin=40e9, fmax=120e9, by_line=1)
17ws.absorption_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
18ws.WignerInit()
19
20# %% Use the automatic agenda setter for propagation matrix calculations
21ws.propagation_matrix_agendaAuto()
22
23# %% Grids and planet
24
25ws.surface_fieldSetPlanetEllipsoid(option="Earth")
26ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0
27ws.atmospheric_fieldRead(
28 toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
29)
30ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37")
31
32# %% Add a sun
33
34ws.sunBlackbody()
35ws.sunsAddSun(suns=ws.suns)
36
37# %% Checks and settings
38
39ws.spectral_radiance_unit = "Tb"
40ws.spectral_radiance_space_agendaSet(option="SunOrCosmicBackground")
41ws.spectral_radiance_surface_agendaSet(option="Blackbody")
42
43# %% Core calculations
44
45pos = [90e3, 0, 0]
46zas = np.linspace(0, 2, 21)
47aas = np.linspace(-180, 180, 21)
48res = np.empty((len(zas), len(aas)))
49for iza in range(len(zas)):
50 for iaa in range(len(aas)):
51 los = [zas[iza], aas[iaa]]
52 ws.ray_pathGeometric(pos=pos, los=los, max_step=1000.0)
53 ws.spectral_radianceClearskyEmission()
54 ws.spectral_radianceApplyUnitFromSpectralRadiance()
55 res[iza, iaa] = ws.spectral_radiance[0][0]
56
57# FIXME: Use some sort of Imager for measurement_vector for the above
58
59r, theta = np.meshgrid(zas, np.rad2deg(aas))
60fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))
61ax.contourf(theta, r, res.T)
62
63assert np.allclose(
64 res[::3, ::7],
65 np.array(
66 [
67 [5.74974671e03, 5.67393003e03, 5.42098021e03],
68 [2.73578415e00, 2.73611538e00, 2.73734690e00],
69 [2.74234763e00, 2.74294365e00, 2.74528926e00],
70 [2.76134660e00, 2.76143226e00, 2.76179796e00],
71 [2.79371400e00, 2.79309756e00, 2.78966340e00],
72 [2.83017769e00, 2.83067328e00, 2.82828924e00],
73 [2.85383040e00, 2.85967097e00, 2.87221221e00],
74 ]
75 ),
76)
Zeeman sun scattering
1import pyarts
2import numpy as np
3import matplotlib.pyplot as plt
4
5ws = pyarts.workspace.Workspace()
6
7# %% Sampled frequency range
8
9ws.frequency_grid = [pyarts.arts.convert.wavelen2freq(700e-9)]
10
11# %% Species and line absorption
12
13ws.absorption_speciesSet(species=["O2-66"])
14ws.ReadCatalogData()
15ws.absorption_bandsSelectFrequency(fmin=40e9, fmax=120e9, by_line=1)
16ws.absorption_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
17ws.WignerInit()
18
19# %% Use the automatic agenda setter for propagation matrix calculations
20ws.propagation_matrix_agendaAuto()
21
22# %% Grids and planet
23
24ws.surface_fieldSetPlanetEllipsoid(option="Earth")
25ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0
26ws.atmospheric_fieldRead(
27 toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
28)
29ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37")
30
31# %% Add a sun
32
33ws.sunBlackbody()
34ws.sunsAddSun(suns=ws.suns)
35
36# %% Checks and settings
37
38ws.spectral_radiance_unit = "Tb"
39ws.spectral_radiance_space_agendaSet(option="SunOrCosmicBackground")
40ws.spectral_radiance_surface_agendaSet(option="Blackbody")
41ws.ray_path_observer_agendaSet(option="Geometric")
42ws.propagation_matrix_scattering_agendaSet(option="AirSimple")
43
44# %% Core calculations
45pos = [90e3, 0, 0]
46zas = np.linspace(0, 5, 21)
47aas = np.linspace(-180, 180, 21)
48res = np.empty((len(zas), len(aas)))
49for iza in range(len(zas)):
50 for iaa in range(len(aas)):
51 los = [zas[iza], aas[iaa]]
52 ws.ray_pathGeometric(pos=pos, los=los, max_step=1000.0)
53 ws.ray_path_suns_pathFromPathObserver(just_hit=1)
54 ws.spectral_radianceClearskyRayleighScattering()
55 ws.spectral_radianceApplyUnitFromSpectralRadiance()
56 res[iza, iaa] = ws.spectral_radiance[0][0]
57
58# FIXME: Use some sort of Imager for measurement_vector for the above
59
60r, theta = np.meshgrid(zas, np.rad2deg(aas))
61fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))
62ax.contourf(theta, r, res.T)
63
64assert np.allclose(
65 res[1::3, 1::7],
66 np.array(
67 [
68 [5771.99999999, 5771.99999999, 5771.99999999],
69 [535.42850284, 535.42853138, 535.42863128],
70 [551.51083023, 551.51086978, 551.51097854],
71 [562.27742756, 562.27747035, 562.27756852],
72 [570.47208475, 570.47212, 570.47223188],
73 [577.13157683, 577.13161927, 577.13172898],
74 [582.7648976, 582.76494322, 582.76505392],
75 ]
76 ),
77)