Zeeman

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)
        angles.append([ang.eta, ang.theta])

fig.tight_layout()

if "ARTS_HEADLESS" not in os.environ:
    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.frequency_grid = np.linspace(-50e6, 50e6, 1001) + line_f0

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

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

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

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

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

# %% Show results
fig, ax = plt.subplots()
ax.plot((ws.frequency_grid - line_f0) / 1e6, ws.spectral_radiance + 0)
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_radiance[::100],
    np.array(
        [[2.27784836e+02,  4.26114598e-04,  1.02751718e-04,  5.69266704e-02],
         [2.30863855e+02,  6.59892211e-04,  1.59299017e-04,  7.04138026e-02],
         [2.34806524e+02,  1.16821510e-03,  2.82508836e-04,  9.33835996e-02],
         [2.40360807e+02,  2.61610579e-03,  6.34821344e-04,  1.40029927e-01],
         [2.49782423e+02,  9.74658804e-03,  2.38967070e-03,  2.69735134e-01],
         [2.09027223e+02,  2.38234847e+01,  1.64634819e+00,  5.44981208e-06],
         [2.49781856e+02,  9.75012917e-03,  2.39056843e-03, -2.69812587e-01],
         [2.40359714e+02,  2.61791138e-03,  6.35266868e-04, -1.40106505e-01],
         [2.34804947e+02,  1.16941576e-03,  2.82802436e-04, -9.34596007e-02],
         [2.30861804e+02,  6.60797190e-04,  1.59519308e-04, -7.04898138e-02],
         [2.27782319e+02,  4.26846028e-04,  1.02929281e-04, -5.70026325e-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()

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

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

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

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

# %% Checks and settings
ws.spectral_radiance_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_vectorFromSensor()

# %% Show results
fig, ax = plt.subplots()
ax.plot((ws.frequency_grid - line_f0) / 1e6, ws.measurement_vector)
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_vector[::100],
    np.array(
        [227.85626444, 230.93431141, 234.89998492, 240.50101182,
         250.05274234, 210.84064948, 249.51259663, 240.21977571,
         234.7115623, 230.79135556, 227.73970923]
    ),
)

(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.frequency_grid = np.linspace(-50e6, 50e6, nf) + line_f0

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

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

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

# %% Checks and settings
ws.spectral_radiance_space_agendaSet(option="Transmission")
ws.spectral_radiance_surface_agendaSet(option="Transmission")

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

# %% Show results
fig, ax = plt.subplots()
ax.semilogy((ws.frequency_grid - line_f0) / 1e6, ws.spectral_radiance)
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_radiance[::100],
    np.array(
        [[3.48824591e-06, -1.47203439e-10, -3.54155511e-11, -4.59223182e-08],
         [1.71707407e-06, -1.11587582e-10, -2.68557096e-11, -2.78282572e-08],
         [6.98942589e-07, -7.92932217e-11, -1.90895244e-11, -1.48129013e-08],
         [2.02345668e-07, -5.04729873e-11, -1.21507664e-11, -6.26440637e-09],
         [2.59661793e-08, -2.51846307e-11, -6.03516645e-12, -1.54349340e-09],
         [6.89347083e-13,  6.83290031e-13, -5.69298697e-14,  7.52268204e-18],
         [2.59895049e-08, -2.52264844e-11, -6.04518733e-12,  1.54841075e-09],
         [2.02661517e-07, -5.06298315e-11, -1.21887251e-11,  6.29847485e-09],
         [7.00441394e-07, -7.96478027e-11, -1.91752558e-11,  1.49217675e-08],
         [1.72166873e-06, -1.12232853e-10, -2.70115093e-11,  2.80800108e-08],
         [3.49929916e-06, -1.48242055e-10, -3.56660238e-11,  4.64085834e-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()

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

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

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

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

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

# %% Checks and settings
ws.spectral_radiance_transform_operatorSet(option="Tb")
ws.spectral_radiance_space_agendaSet(option="SunOrCosmicBackground")
ws.spectral_radiance_surface_agendaSet(option="Blackbody")

# %% Core calculations
pos = [90e3, 0, 0]
zas = np.linspace(0, 2, 21)
aas = np.linspace(-180, 180, 21)
res = np.empty((len(zas), len(aas)))
for iza in range(len(zas)):
    for iaa in range(len(aas)):
        los = [zas[iza], aas[iaa]]
        ws.ray_pathGeometric(pos=pos, los=los, max_step=1000.0)
        ws.spectral_radianceClearskyEmission()
        ws.spectral_radianceApplyUnitFromSpectralRadiance()
        res[iza, iaa] = ws.spectral_radiance[0][0]

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

r, theta = np.meshgrid(zas, np.rad2deg(aas))
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(
        [[5333.26510189, 5333.26510189, 5333.26510189],
         [  17.54346643,   17.42489325,   17.45953778],
         [  17.61367450,   17.37693775,   17.44655483],
         [  17.68367603,   17.32918821,   17.43410485],
         [  17.75346986,   17.28164635,   17.42218858],
         [  17.82305482,   17.23431393,   17.41080683],
         [  17.89242977,   17.18719274,   17.39996041]]
    ),
)

(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.frequency_grid = [pyarts.arts.convert.wavelen2freq(700e-9)]

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

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

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

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

# %% Checks and settings
ws.spectral_radiance_transform_operator = "Tb"
ws.spectral_radiance_space_agendaSet(option="SunOrCosmicBackground")
ws.spectral_radiance_surface_agendaSet(option="Blackbody")
ws.ray_path_observer_agendaSetGeometric()
ws.propagation_matrix_scattering_agendaSet(option="AirSimple")

# %% Core calculations
pos = [90e3, 0, 0]
zas = np.linspace(0, 5, 21)
aas = np.linspace(-180, 180, 21)
res = np.empty((len(zas), len(aas)))
for iza in range(len(zas)):
    for iaa in range(len(aas)):
        los = [zas[iza], aas[iaa]]
        ws.ray_pathGeometric(pos=pos, los=los, max_step=1000.0)
        ws.ray_path_suns_pathFromPathObserver(just_hit=1)
        ws.spectral_radianceClearskyRayleighScattering()
        ws.spectral_radianceApplyUnitFromSpectralRadiance()
        res[iza, iaa] = ws.spectral_radiance[0][0]

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

r, theta = np.meshgrid(zas, np.rad2deg(aas))
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.99991010, 5771.99991010, 5771.99991010],
              [ 643.69008277,  643.69008277,  643.69008277],
              [ 643.69030757,  643.69030764,  643.69030764],
              [ 643.69065708,  643.69065715,  643.69065728],
              [ 643.69113428,  643.69113428,  643.69113441],
              [ 643.69174257,  643.69174270,  643.69174290],
              [ 643.69248695,  643.69248708,  643.69248747]]),
)

(Source code, svg, pdf)

_images/6-zeeman-sun-scattering.svg