Zeeman

Zeeman radiative-transfer examples

These examples exercise polarized clear-sky radiative transfer around the 118.75 GHz oxygen line. In addition to geometry, sensors, refraction, solar backgrounds, and adaptive paths, the directory compares the available within-layer propagation and source approximations selected by ws.rte_option.

The relevant examples are:

  • 2-zeeman.py uses constant;

  • 9-zeeman-linear-src.py uses lintau;

  • 10-zeeman-linear-prop.py uses linprop;

  • 11-zeeman-magnus.py uses magop; and

  • 12-zeeman-magnus-linear-src.py uses magop_linsrc.

magop represents the propagation matrix as linear across each layer and includes the first commutator correction of the Magnus expansion, while using the endpoint-average source. magop_linsrc additionally represents the source as linear across the layer using the augmented Magnus source operator. Both options support polarized analytical Jacobians.

The regression arrays sampled in these examples are intentional: the Magnus results need not equal constant, lintau, or linprop when adjacent polarized propagation matrices do not commute. For a constant or commuting propagation matrix, magop reduces to the ordinary endpoint-average matrix exponential.

Zeeman geometry

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# %% Get magnetic field at a sample position
ell = pyarts.arts.planets.Earth.ellipsoid
pos = [50e3, 0, 0]
mag = pyarts.arts.igrf(pos, ell, t="2000-03-11 14:39:37")

# %% Setup figure for multiple subplots
N = 3
M = 4
fig = plt.figure(figsize=(M * 8, N * 8))

# %% Store computed angles for later verification
angles = []

# %% loop over different LOS directions
for i in range(N):
    for j in range(M):
        los = [np.linspace(0, 180, N)[i], np.linspace(0, 360, M)[j]]

        # Setup 3D subplot
        ax = fig.add_subplot(N, M, i*M + j + 1, projection='3d')

        # Plot and store angles
        ang = pyarts.arts.zeeman.MagneticAngles(mag, los)
        pyarts.plots.MagneticAngles.plot(ang, fig=fig, ax=ax, N=50)
        ax.set_ylim(-1.2, 1.2)
        ax.set_xlim(-1.2, 1.2)
        ax.set_zlim(-1.2, 1.2)
        angles.append([ang.eta, ang.theta])

if "ARTS_HEADLESS" not in os.environ:
    plt.tight_layout()
    plt.show()

assert np.allclose(angles,
                   np.array([[3.02347622,  1.07933684],
                             [0.92908112,  1.07933684],
                             [-1.16531399,  1.07933684],
                             [3.02347622,  1.07933684],
                             [0.21669967,  0.50432248],
                             [0.98174228,  2.12671818],
                             [-1.04334554,  1.92599378],
                             [0.21669967,  0.50432248],
                             [0.11811643,  2.06225581],
                             [2.21251154,  2.06225581],
                             [-1.97627867,  2.06225581],
                             [0.11811643,  2.06225581]])), "Angles do not match expected values."

(Source code, svg, pdf)

_images/1-zeeman-geometry.svg

Zeeman

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()

# %% Sampled frequency range
line_f0 = 118750348044.712
ws.freq_grid = np.linspace(-50e6, 50e6, 1001) + line_f0

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldSchmidthFieldFromIGRF(time="2000-03-11 14:39:37")

# %% Checks and settings
ws.spectral_rad_transform_operatorSet(option="Tb")

# %% Core calculations
pos = [100e3, 0, 0]
los = [180.0, 0.0]
ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
ws.rte_option = "constant"
ws.spectral_radClearskyEmission()
ws.spectral_radApplyUnitFromSpectralRadiance()

# %% Show results
fig, ax = pyarts.plot(ws.spectral_rad, freqs=(
    ws.freq_grid - line_f0) / 1e6)
[a.set_xlabel("Frequency offset [MHz]") for a in ax.flatten()]
[a.set_ylabel("Spectral radiance [K]") for a in ax.flatten()]
fig.suptitle(f"Zeeman effect of {round(line_f0 / 1e6)} MHz O$_2$ line")

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

# %% Test

assert np.allclose(
    ws.spectral_rad[::100],
    np.array(
        [[ 2.27651282e+02,  4.26102060e-04,  1.02747352e-04,  5.68792738e-02],
         [ 2.30728866e+02,  6.60200882e-04,  1.59372185e-04,  7.04074714e-02],
         [ 2.34671712e+02,  1.16870303e-03,  2.82624380e-04,  9.34017354e-02],
         [ 2.40226524e+02,  2.61619467e-03,  6.34833158e-04,  1.40041532e-01],
         [ 2.49649981e+02,  9.75455461e-03,  2.39140470e-03,  2.69997141e-01],
         [ 2.09901775e+02,  2.41594007e+01,  1.73753634e+00,  5.52462292e-06],
         [ 2.49649415e+02,  9.75809646e-03,  2.39230231e-03, -2.70074598e-01],
         [ 2.40225433e+02,  2.61800016e-03,  6.35278628e-04, -1.40118090e-01],
         [ 2.34670137e+02,  1.16990381e-03,  2.82918000e-04, -9.34777147e-02],
         [ 2.30726819e+02,  6.61105602e-04,  1.59592410e-04, -7.04834216e-02],
         [ 2.27648772e+02,  4.26832612e-04,  1.02924697e-04, -5.69551084e-02]]
    ),
), "Values have drifted from expected results in spectral radiance"

(Source code, svg, pdf)

_images/2-zeeman.svg

Zeeman sensor

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()
ws.rte_option = "constant"

# %% Sampled frequency range
line_f0 = 118750348044.712
ws.freq_grid = np.linspace(-50e6, 50e6, 1001) + line_f0

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldIGRF(time="2000-03-11 14:39:37")

# %% Checks and settings
ws.spectral_rad_transform_operator = "Tb"
ws.ray_path_observer_agendaSetGeometric()

# %% Set up a sensor with Gaussian standard deviation channel widths on individual frequency ranges
pos = [100e3, 0, 0]
los = [180.0, 0.0]
ws.measurement_sensorSimpleGaussian(std=1e5, pos=pos, los=los, pol="RC")

# %% Core calculations
ws.measurement_vecFromSensor()

# %% Show results
fig, ax = pyarts.plot(ws.measurement_vec, xgrid=(
    ws.freq_grid - line_f0) / 1e6)
ax.set_xlabel("Frequency offset [MHz]")
ax.set_ylabel("Spectral radiance [K]")
ax.set_title(
    f"Zeeman effect of {round(line_f0 / 1e6)} MHz O$_2$ line with Gaussian channels on individual grids"
)

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

# %% Test
assert np.allclose(
    ws.measurement_vec[::100],
    np.array(
        [227.72265120, 230.79931575, 234.76519053, 240.36674038,
         249.92056686, 211.81135165, 249.37989662, 240.08548277,
         234.57673412, 230.65637644, 227.60619667]
    ),
)

(Source code, svg, pdf)

_images/3-zeeman-sensor.svg

Zeeman transmission

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()

# %% Sampled frequency range
line_f0 = 118750348044.712
nf = 1001
ws.freq_grid = np.linspace(-50e6, 50e6, nf) + line_f0

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldIGRF(time="2000-03-11 14:39:37")

# %% Checks and settings
ws.spectral_rad_space_agendaSet(option="Transmission")
ws.spectral_rad_surface_agendaSet(option="Transmission")

# %% Core calculations
pos = [100e3, 0, 0]
los = [180.0, 0.0]
ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
ws.spectral_radClearskyTransmission()

# %% Show results
fig, ax = pyarts.plot(ws.spectral_rad, freqs=(
    ws.freq_grid - line_f0) / 1e6, component='I')
ax.set_yscale('log')
ax.set_xlabel("Frequency offset [MHz]")
ax.set_ylabel("Spectral radiance [K]")
ax.set_title(f"Zeeman effect of {round(line_f0 / 1e6)} MHz O$_2$ line")

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

# %% Test
assert np.allclose(
    ws.spectral_rad[::100],
    np.array(
        [[ 3.55416974e-06, -1.49626896e-10, -3.59945621e-11, -4.65675860e-08],
         [ 1.75504614e-06, -1.14112664e-10, -2.74625108e-11, -2.83292389e-08],
         [ 7.16577429e-07, -8.14996162e-11, -1.96228377e-11, -1.51461032e-08],
         [ 2.07916647e-07, -5.19869800e-11, -1.25199070e-11, -6.43052931e-09],
         [ 2.67056326e-08, -2.58535066e-11, -6.20359298e-12, -1.58692924e-09],
         [ 7.87383592e-13,  7.80210615e-13, -6.46005354e-14,  8.84600141e-18],
         [ 2.67295620e-08, -2.58964360e-11, -6.21391407e-12,  1.59197636e-09],
         [ 2.08240338e-07, -5.21480127e-11, -1.25589127e-11,  6.46545336e-09],
         [ 7.18110068e-07, -8.18625471e-11, -1.97105925e-11,  1.52574035e-08],
         [ 1.75973088e-06, -1.14770061e-10, -2.76211969e-11,  2.85857722e-08],
         [ 3.56540626e-06, -1.50680035e-10, -3.62484502e-11,  4.70615011e-08]]
    ),
), "Values have drifted from expected results in spectral radiance"

(Source code, svg, pdf)

_images/4-zeeman-transmission.svg

Zeeman sun

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()
ws.rte_option = "constant"

# %% Sampled frequency range
line_f0 = 118750348044.712
ws.freq_grid = [line_f0]

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldIGRF(time="2000-03-11 14:39:37")

# %% Add a sun
ws.sunBlackbody()
ws.suns = [ws.sun]

# %% Checks and settings
ws.spectral_rad_transform_operatorSet(option="Tb")
ws.spectral_rad_space_agendaSet(option="SunOrCosmicBackground")
ws.spectral_rad_surface_agendaSet(option="Blackbody")

# %% Core calculations
pos = [90e3, 0, 0]
zens = np.linspace(0, 2, 21)
azis = np.linspace(-180, 180, 21)
res = np.empty((len(zens), len(azis)))
for izen in range(len(zens)):
    for iazi in range(len(azis)):
        los = [zens[izen], azis[iazi]]
        ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
        ws.spectral_radClearskyEmission()
        ws.spectral_radApplyUnitFromSpectralRadiance()
        res[izen, iazi] = ws.spectral_rad[0][0]

# FIXME: Use some sort of Imager for measurement_vec for the above

r, theta = np.meshgrid(zens, np.rad2deg(azis))
fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))
ax.contourf(theta, r, res.T)

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

assert np.allclose(
    res[::3, ::7],
    np.array(
        [[5356.88051075, 5356.88051075, 5356.88051075],
         [  16.77380467,   16.66091577,   16.69389844],
         [  16.84063115,   16.61524188,   16.68151802],
         [  16.90725049,   16.56975243,   16.6696318 ],
         [  16.9736614 ,   16.52444918,   16.65824055],
         [  17.03986264,   16.47933394,   16.64734506],
         [  17.10585295,   16.43440858,   16.63694619]])
    )

(Source code, svg, pdf)

_images/5-zeeman-sun.svg

Zeeman sun scattering

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()

# %% Sampled frequency range
ws.freq_grid = [pyarts.arts.convert.wavelen2freq(700e-9)]

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldIGRF(time="2000-03-11 14:39:37")

# %% Add a sun
ws.sunBlackbody()
ws.suns = [ws.sun]

# %% Checks and settings
ws.spectral_rad_transform_operator = "Tb"
ws.spectral_rad_space_agendaSet(option="SunOrCosmicBackground")
ws.spectral_rad_surface_agendaSet(option="Blackbody")
ws.ray_path_observer_agendaSetGeometric()
ws.spectral_propmat_scat_agendaSet(option="AirSimple")

# %% Core calculations
pos = [90e3, 0, 0]
zens = np.linspace(0, 5, 21)
azis = np.linspace(-180, 180, 21)
res = np.empty((len(zens), len(azis)))
for izen in range(len(zens)):
    for iazi in range(len(azis)):
        los = [zens[izen], azis[iazi]]
        ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
        ws.ray_path_suns_pathFromPathObserver(just_hit=1)
        ws.spectral_radClearskyRayleighScattering()
        ws.spectral_radApplyUnitFromSpectralRadiance()
        res[izen, iazi] = ws.spectral_rad[0][0]

# FIXME: Use some sort of Imager for measurement_vec for the above

r, theta = np.meshgrid(zens, np.rad2deg(azis))
fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))
ax.contourf(theta, r, res.T)

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

assert np.allclose(
    res[1::3, 1::7],
    np.array([[5771.9999155 , 5771.9999155 , 5771.9999155 ],
              [ 642.44289971,  642.44289971,  642.44289971],
              [ 642.44286517,  642.44286517,  642.44286524],
              [ 642.4428158 ,  642.4428158 ,  642.44281587],
              [ 642.44275627,  642.44275634,  642.44275648],
              [ 642.44269338,  642.44269352,  642.44269372],
              [ 642.44263584,  642.44263598,  642.44263625]]),
)

(Source code, svg, pdf)

_images/6-zeeman-sun-scattering.svg

Zeeman refractive

import os

import pyarts3 as pyarts
import numpy as np
import matplotlib.pyplot as plt

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()

# %% Sampled frequency range

line_f0 = 53.0669e9
ws.freq_grid = np.linspace(-15e6, 15e6, 51) + line_f0

# %% Species and line absorption

ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=52e9, fmax=54e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet

ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=120e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldIGRF()

# %% Checks and settings

ws.spectral_rad_transform_operatorSet(option="Tb")


@pyarts.arts_agenda(ws=ws, fix=True)
def single_propmat_agenda(ws):
    ws.single_propmatInit()
    ws.single_propmatAddVoigtLTE()


# %% Calculate and compare refractive and geometric ray paths
pos = [3571, 46, 7]
los = [20, 90]

res = []
ws.obs_pos = pos
ws.obs_los = los
ws.max_stepsize = 20000.0

# %% Show results
ws.ray_point_back_propagation_agendaSet(option="GeometricStepwise")
ws.spectral_radClearskyEmissionFrequencyDependentPropagation(max_tau=1e-2)
ws.spectral_radApplyUnitFromSpectralRadiance(ray_path=ws.spectral_ray_path[0])
geometric = ws.spectral_rad * 1.0
ws.ray_point_back_propagation_agendaSet(option="RefractiveStepwise")
ws.spectral_radClearskyEmissionFrequencyDependentPropagation(max_tau=1e-2)
ws.spectral_radApplyUnitFromSpectralRadiance(ray_path=ws.spectral_ray_path[0])
refractive = ws.spectral_rad * 1.0

# %% Show results

if "ARTS_HEADLESS" not in os.environ:
    fig, ax = plt.subplots(2, 2, figsize=(10, 8))
    freqs = ws.freq_grid / 1e9
    pyarts.plots.StokvecVector.plot(
        geometric - refractive, fig=fig, ax=ax, freqs=freqs)
    for a in ax.flatten():
        a.set_xlabel("Frequency offset [MHz]")
        a.set_ylabel("Spectral radiance [K]")
    ax[0, 0].set_title("Stokes I")
    ax[0, 1].set_title("Stokes Q")
    ax[1, 0].set_title("Stokes U")
    ax[1, 1].set_title("Stokes V")
    fig.suptitle(
        f"Difference between geometric and refractive paths around the {round(line_f0 / 1e6)} MHz O$_2$ line")
    fig.tight_layout()
    plt.show()

# %% Tests

assert np.allclose(geometric.flatten()[::21], [1.25486103e+02, -5.56777144e-04, -2.63816009e-03,  6.12431001e-01,
                                               1.38429184e+02,  1.04834522e-01, -2.21549573e-02, -5.18205673e-01,
                                               1.29365186e+02, -4.48679095e-04])

assert np.allclose(refractive.flatten()[::21], [1.25552512e+02, -5.61580828e-04, -2.61219720e-03,  6.11932940e-01,
                                                1.38488011e+02,  1.04448633e-01, -2.18946938e-02, -5.19985291e-01,
                                                1.29439037e+02, -4.53301859e-04])

assert not np.allclose(geometric, refractive)

(Source code)

Zeeman adaptive path

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts
from copy import deepcopy as copy

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()

# %% Sampled frequency range
line_f0 = 118750348044.712
ws.freq_grid = np.linspace(-50e6, 50e6, 1001) + line_f0

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldSchmidthFieldFromIGRF(time="2000-03-11 14:39:37")

# %% Checks and settings
ws.spectral_rad_transform_operatorSet(option="Tb")

# %% Core calculations
pos = [100e3, 0, 0]
los = [180.0, 0.0]
ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
ws.ray_pointBackground()
ws.spectral_rad_bkgAgendasAtEndOfPath()
ws.atm_pathFromPath()
ws.freq_grid_pathFromPath()
ws.spectral_propmat_pathFromPath()
ws.spectral_radSetToBackground()
ws.spectral_radSinglePathEmissionFrequencyLoop()
ws.spectral_radApplyUnitFromSpectralRadiance()

srad0 = copy(ws.spectral_rad)
path0 = copy(ws.ray_path)

ws.spectral_propmat_pathAddAdaptiveHalfPath(
    max_stepsize=100., max_tau=0.05, cutoff_tau=3.0)
ws.spectral_radSetToBackground()
ws.spectral_radSinglePathEmissionFrequencyLoop()
ws.spectral_radApplyUnitFromSpectralRadiance()

freqs = (ws.freq_grid - line_f0) / 1e6
f, a = pyarts.plot(srad0, freqs=freqs, label="Regular path")
[a.set_xlabel("Frequency offset [MHz]") for a in a.flatten()]
[a.set_ylabel("Spectral radiance [K]") for a in a.flatten()]
f.suptitle(f"Zeeman effect of {round(line_f0 / 1e6)} MHz O$_2$ line")
pyarts.plot(ws.spectral_rad, freqs=freqs, fig=f, ax=a, label="Adaptive path")
[a.legend() for a in a.flatten()]

fig, ax = plt.subplots(1, 1)
d0 = path0.distances(ws.surf_field.ellipsoid)
a0 = [p.pos[0] / 1e3 for p in path0]
ax.plot(d0, a0[:-1], label="Regular path")
d1 = ws.ray_path.distances(ws.surf_field.ellipsoid)
a1 = [p.pos[0] / 1e3 for p in ws.ray_path]
ax.plot(d1, a1[:-1], label="Adaptive path")
ax.legend()
ax.set_xlabel("Distance along path [m]")
ax.set_ylabel("Altitude [km]")
fig.suptitle("Ray paths")

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

assert len(ws.ray_path) > len(path0)  # More points in adaptive path
assert np.allclose([np.sum(d0)], [np.sum(d1)])
assert np.allclose(sorted(a0, reverse=True), a0)
assert np.allclose(sorted(a1, reverse=True), a1)

(Source code)

_images/8-zeeman-adaptive-path_00.svg

Fig. 1 (svg, pdf)

_images/8-zeeman-adaptive-path_01.svg

Fig. 2 (svg, pdf)

Zeeman linear src

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()

# %% Sampled frequency range
line_f0 = 118750348044.712
ws.freq_grid = np.linspace(-50e6, 50e6, 1001) + line_f0

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldSchmidthFieldFromIGRF(time="2000-03-11 14:39:37")

# %% Checks and settings
ws.spectral_rad_transform_operatorSet(option="Tb")

# %% Core calculations
pos = [100e3, 0, 0]
los = [180.0, 0.0]
ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
ws.rte_option = "lintau"
ws.spectral_radClearskyEmission()
ws.spectral_radApplyUnitFromSpectralRadiance()

# %% Show results
fig, ax = pyarts.plot(ws.spectral_rad, freqs=(
    ws.freq_grid - line_f0) / 1e6)
[a.set_xlabel("Frequency offset [MHz]") for a in ax.flatten()]
[a.set_ylabel("Spectral radiance [K]") for a in ax.flatten()]
fig.suptitle(f"Zeeman effect of {round(line_f0 / 1e6)} MHz O$_2$ line")

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

# %% Test

assert np.allclose(
    ws.spectral_rad[::100],
    np.array(
        [[ 2.27693776e+02,  4.25360914e-04,  1.02565807e-04,  5.68301063e-02],
         [ 2.30768819e+02,  6.59206165e-04,  1.59126872e-04,  7.03542409e-02],
         [ 2.34709073e+02,  1.16721685e-03,  2.82253283e-04,  9.33505045e-02],
         [ 2.40261448e+02,  2.61334833e-03,  6.34105014e-04,  1.39994065e-01],
         [ 2.49682223e+02,  9.74387048e-03,  2.38851618e-03,  2.69924023e-01],
         [ 2.09873451e+02,  2.41579480e+01,  1.73827372e+00,  5.54548559e-06],
         [ 2.49681656e+02,  9.74741170e-03,  2.38941342e-03, -2.70001562e-01],
         [ 2.40260356e+02,  2.61515328e-03,  6.34550292e-04, -1.40070694e-01],
         [ 2.34707497e+02,  1.16841699e-03,  2.82546721e-04, -9.34265488e-02],
         [ 2.30766770e+02,  6.60110231e-04,  1.59346921e-04, -7.04302667e-02],
         [ 2.27691263e+02,  4.26090893e-04,  1.02743001e-04, -5.69060490e-02]]
    ),
), "Values have drifted from expected results in spectral radiance"

(Source code, svg, pdf)

_images/9-zeeman-linear-src.svg

Zeeman linear prop

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()

# %% Sampled frequency range
line_f0 = 118750348044.712
ws.freq_grid = np.linspace(-50e6, 50e6, 1001) + line_f0

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldSchmidthFieldFromIGRF(time="2000-03-11 14:39:37")

# %% Checks and settings
ws.spectral_rad_transform_operatorSet(option="Tb")

# %% Core calculations
pos = [100e3, 0, 0]
los = [180.0, 0.0]
ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
ws.rte_option = "linprop"
ws.spectral_radClearskyEmission()
ws.spectral_radApplyUnitFromSpectralRadiance()

# %% Show results
fig, ax = pyarts.plot(ws.spectral_rad, freqs=(
    ws.freq_grid - line_f0) / 1e6)
[a.set_xlabel("Frequency offset [MHz]") for a in ax.flatten()]
[a.set_ylabel("Spectral radiance [K]") for a in ax.flatten()]
fig.suptitle(f"Zeeman effect of {round(line_f0 / 1e6)} MHz O$_2$ line")

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

# %% Test

assert np.allclose(
    ws.spectral_rad[::100],
    np.array(
        [[ 2.27655097e+02,  4.25941643e-04,  1.02708595e-04,  5.68661397e-02],
         [ 2.30731823e+02,  6.60022963e-04,  1.59329204e-04,  7.03869074e-02],
         [ 2.34673691e+02,  1.16858889e-03,  2.82596751e-04,  9.33841959e-02],
         [ 2.40228067e+02,  2.61655084e-03,  6.34919268e-04,  1.40046082e-01],
         [ 2.49653926e+02,  9.76155919e-03,  2.39312012e-03,  2.70138698e-01],
         [ 2.09894520e+02,  2.41722890e+01,  1.73860327e+00,  5.52740752e-06],
         [ 2.49653359e+02,  9.76510257e-03,  2.39401813e-03, -2.70216174e-01],
         [ 2.40226976e+02,  2.61835630e-03,  6.35364732e-04, -1.40122633e-01],
         [ 2.34672116e+02,  1.16978940e-03,  2.82890308e-04, -9.34601655e-02],
         [ 2.30729776e+02,  6.60927425e-04,  1.59549365e-04, -7.04628681e-02],
         [ 2.27652586e+02,  4.26672085e-04,  1.02885912e-04, -5.69420287e-02]]
    ),
), "Values have drifted from expected results in spectral radiance"

(Source code, svg, pdf)

_images/10-zeeman-linear-prop.svg

Zeeman magnus

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()

# %% Sampled frequency range
line_f0 = 118750348044.712
ws.freq_grid = np.linspace(-50e6, 50e6, 1001) + line_f0

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldSchmidthFieldFromIGRF(time="2000-03-11 14:39:37")

# %% Checks and settings
ws.spectral_rad_transform_operatorSet(option="Tb")

# %% Core calculations
pos = [100e3, 0, 0]
los = [180.0, 0.0]
ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
ws.rte_option = "magop"
ws.spectral_radClearskyEmission()
ws.spectral_radApplyUnitFromSpectralRadiance()

# %% Show results
fig, ax = pyarts.plot(ws.spectral_rad, freqs=(ws.freq_grid - line_f0) / 1e6)
[a.set_xlabel("Frequency offset [MHz]") for a in ax.flatten()]
[a.set_ylabel("Spectral radiance [K]") for a in ax.flatten()]
fig.suptitle(f"Zeeman effect of {round(line_f0 / 1e6)} MHz O$_2$ line")

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

# %% Test

assert np.allclose(
    ws.spectral_rad[::100],
    np.array(
        [
            [2.276512821092e02, 4.261021957888e-04, 1.027468021420e-04, 5.687927378416e-02],
            [2.307288660365e02, 6.602011364976e-04, 1.593711565846e-04, 7.040747139975e-02],
            [2.346717117551e02, 1.168703614880e-03, 2.826220439971e-04, 9.340173537775e-02],
            [2.402265237320e02, 2.616196566095e-03, 6.348256699340e-04, 1.400415316211e-01],
            [2.496499814427e02, 9.754569416941e-03, 2.391349108080e-03, 2.699971405358e-01],
            [2.099028163397e02, 2.416037218539e01, 1.740559369663e00, 5.523843341848e-06],
            [2.496494147479e02, 9.758111279439e-03, 2.392246667185e-03, -2.700745982318e-01],
            [2.402254325160e02, 2.618002063414e-03, 6.352711243096e-04, -1.401180900068e-01],
            [2.346701369088e02, 1.169904391432e-03, 2.829156569959e-04, -9.347771470715e-02],
            [2.307268185042e02, 6.611058571622e-04, 1.595913768142e-04, -7.048342160883e-02],
            [2.276487716773e02, 4.268327476638e-04, 1.029241438602e-04, -5.695510840674e-02],
        ]
    ),
), "Magnus Zeeman spectral radiance has drifted from the expected values"

(Source code, svg, pdf)

_images/11-zeeman-magnus.svg

Zeeman magnus linear src

import os

import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts

# Download catalogs
pyarts.data.download()

ws = pyarts.workspace.Workspace()

# %% Sampled frequency range
line_f0 = 118750348044.712
ws.freq_grid = np.linspace(-50e6, 50e6, 1001) + line_f0

# %% Species and line absorption
ws.abs_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
ws.WignerInit()

# %% Use the automatic agenda setter for propagation matrix calculations
ws.spectral_propmat_agendaAuto()

# %% Grids and planet
ws.surf_fieldPlanet(option="Earth")
ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atm_fieldRead(
    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atm_fieldSchmidthFieldFromIGRF(time="2000-03-11 14:39:37")

# %% Checks and settings
ws.spectral_rad_transform_operatorSet(option="Tb")

# %% Core calculations
pos = [100e3, 0, 0]
los = [180.0, 0.0]
ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
ws.rte_option = "magop_linsrc"
ws.spectral_radClearskyEmission()
ws.spectral_radApplyUnitFromSpectralRadiance()

# %% Show results
fig, ax = pyarts.plot(ws.spectral_rad, freqs=(ws.freq_grid - line_f0) / 1e6)
[a.set_xlabel("Frequency offset [MHz]") for a in ax.flatten()]
[a.set_ylabel("Spectral radiance [K]") for a in ax.flatten()]
fig.suptitle(f"Zeeman effect of {round(line_f0 / 1e6)} MHz O$_2$ line")

if "ARTS_HEADLESS" not in os.environ:
    plt.show()

# %% Test

assert np.allclose(
    ws.spectral_rad[::100],
    np.array(
        [
            [2.276551235399e02, 4.259411802252e-04, 1.027078921965e-04, 5.686609432459e-02],
            [2.307318484287e02, 6.600222918109e-04, 1.593279395564e-04, 7.038684472572e-02],
            [2.346737153328e02, 1.168587845157e-03, 2.825940092634e-04, 9.338410284805e-02],
            [2.402280899194e02, 2.616549259585e-03, 6.349109213630e-04, 1.400459368601e-01],
            [2.496539453642e02, 9.761561980270e-03, 2.393061571382e-03, 2.701384329393e-01],
            [2.098955622505e02, 2.417326246041e01, 1.741633379836e00, 5.526631596808e-06],
            [2.496533784517e02, 9.765105367691e-03, 2.393959518344e-03, -2.702159091422e-01],
            [2.402269987255e02, 2.618354729037e-03, 6.353563693580e-04, -1.401224878107e-01],
            [2.346721406968e02, 1.169788358766e-03, 2.828875578587e-04, -9.346007246592e-02],
            [2.307298009667e02, 6.609267549322e-04, 1.595480962493e-04, -7.046280541539e-02],
            [2.276526123549e02, 4.266716225771e-04, 1.028852062035e-04, -5.694198344678e-02],
        ]
    ),
), "Linear-source Magnus Zeeman spectral radiance has drifted from the expected values"

(Source code, svg, pdf)

_images/12-zeeman-magnus-linear-src.svg