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)