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
)
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
)
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
)
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
)
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
)
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
)