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.27786279e02, -2.23644583e-04, -3.77015913e-04, 5.69269522e-02],
60 [2.30865312e02, -3.48096654e-04, -5.82873607e-04, 7.04140485e-02],
61 [2.34807995e02, -6.21529170e-04, -1.02888198e-03, 9.33839855e-02],
62 [2.40362302e02, -1.41676334e-03, -2.28982732e-03, 1.40030558e-01],
63 [2.49783915e02, -5.60528357e-03, -8.33754971e-03, 2.69733178e-01],
64 [2.07245890e02, -4.30794421e00, -2.14467490e01, 1.24566607e-05],
65 [2.49783382e02, -5.60754703e-03, -8.34040843e-03, -2.69809821e-01],
66 [2.40361281e02, -1.41779460e-03, -2.29136074e-03, -1.40106384e-01],
67 [2.34806526e02, -6.22190182e-04, -1.02991615e-03, -9.34592546e-02],
68 [2.30863403e02, -3.48585873e-04, -5.83658317e-04, -7.04893428e-02],
69 [2.27783939e02, -2.24035725e-04, -3.77652642e-04, -5.70022211e-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_sensorGaussian(
45 f0_fwhm_df=[[f, 1e5, 1e4] for f in ws.frequency_grid],
46 pos=pos,
47 los=los,
48 pol="RC",
49)
50
51# %% Core calculations
52
53result = pyarts.arts.Vector()
54result_jac = pyarts.arts.Matrix()
55ws.measurement_vectorFromSensor(result, result_jac)
56
57# %% Show results
58
59plt.plot((ws.frequency_grid - line_f0) / 1e6, result)
60plt.xlabel("Frequency offset [MHz]")
61plt.ylabel("Spectral radiance [K]")
62plt.title(
63 f"Zeeman effect of {round(line_f0/1e6)} MHz O$_2$ line with Gaussian channels on individual grids"
64)
65
66
67# %% Test
68
69assert np.allclose(
70 result[::100],
71 np.array(
72 [
73 227.84323245,
74 230.93576779,
75 234.90145621,
76 240.50250701,
77 250.05423085,
78 209.93307663,
79 249.51412418,
80 240.22134231,
81 234.71314135,
82 230.79295502,
83 227.72696315,
84 ]
85 ),
86)
87
88
89# %% Set up a sensor with a fixed grid size
90
91pos = [100e3, 0, 0]
92los = [180.0, 0.0]
93ws.measurement_sensorGaussianFrequencyGrid(
94 f0_fwhm=[[f, 1e5] for f in ws.frequency_grid],
95 pos=pos,
96 los=los,
97 pol="RC",
98)
99
100# %% Core calculations
101
102result = pyarts.arts.Vector()
103result_jac = pyarts.arts.Matrix()
104ws.measurement_vectorFromSensor(result, result_jac)
105
106# %% Show results
107
108plt.plot((ws.frequency_grid - line_f0) / 1e6, result)
109plt.xlabel("Frequency offset [MHz]")
110plt.ylabel("Spectral radiance [K]")
111plt.title(
112 f"Zeeman effect of {round(line_f0/1e6)} MHz O$_2$ line with Gaussian channels on a single grid"
113)
114
115
116# %% Test
117
118assert np.allclose(
119 result[::100],
120 np.array(
121 [
122 227.85768946,
123 230.9357677,
124 234.90145604,
125 240.50250662,
126 250.05422954,
127 209.91116972,
128 249.51412294,
129 240.22134193,
130 234.71314118,
131 230.79295493,
132 227.74131132,
133 ]
134 ),
135)
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)