Absorption cross sections
Lines
"""
This file will showcase how you can load line data from arts-cat-data into the
workspace and do the required setups to perform a simple forward calculations
using this data.
"""
import os
import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts
# Initialize ARTS workspace
ws = pyarts.workspace.Workspace()
"""
The workspace contains an absorption species helper list
that allows it to understand what data and methods you
will want to call.
This example sets the absorption species to all oxygen
and all water isotopologues (in their line-by-line mode).
"""
ws.absorption_speciesSet(species=["O2-66", "H2O-161"])
"""
We now need to load the data.
ARTS comes with a lot of line-by-line data available
in its arts-cat-data repository. The first line below
downloads that data to a local cache and allows you
to run this script.
You can also provide line data in other ways, e.g., by loading HITRAN
or some other format. What is important is to populate the absorption
bands with appropriate data.
"""
pyarts.data.download()
ws.absorption_bandsReadSpeciesSplitCatalog(basename="lines/")
"""
Compute absorption
Now we can use any of the methods available to compute line-by-line absorption.
Below we simply set a simple atmosphere and plot the resulting absorption.
Please see other examples for more details on how to use line-by-line
calculations in increasingly complex ways to solve your problem
"""
# Initialize an atmospheric point
atm = pyarts.arts.AtmPoint()
atm.temperature = 295 # At room temperature
atm.pressure = 1e5 # At 1 bar
atm["O2"] = 0.21 # At 21% atmospheric Oxygen
atm["H2O"] = 0.001 # At 0.1% atmospheric Water Vapor
# Set a frequency range and remove lines outside (to speed up calculations)
freqs = np.linspace(1e9, 1000e9, 1001)
ws.absorption_bands.keep_frequencies(fmax=freqs[-1])
# Use available plotting routines to draw the results
fig, ax = pyarts.plot(
ws.absorption_bands,
freqs=freqs,
atm=atm,
mode='important fill isotopes',
min_pm=1e-7)
ax.legend(bbox_to_anchor=(1.05, 1))
ax.set_xlim(freqs[0], freqs[-1])
ax.set_ylim(ymin=1e-7, ymax=1)
ax.set_xticks(np.linspace(0, 1000e9, 11), np.linspace(0, 1000, 11))
ax.set_yscale('log')
ax.set_xlabel("Frequency [GHz]")
ax.set_ylabel("Absorption [1/m]")
ax.set_title("Absorption by species")
if "ARTS_HEADLESS" not in os.environ:
plt.show()
(Source code, svg, pdf)
Cia
"""
This file will showcase how you can load collision-induced absorption (CIA)
data from arts-cat-data into the workspace and do the required setups to
perform a simple forward calculations using this data
Note that this example presumes that you have set the environment variable
ARTS_DATA_PATH to contain a path to a local copy of both arts-cat-data and
arts-xml-data before you import pyarts3. Please check that this is the case
if the example does not work for you. You can easily check if this path is
set by adding the following two lines at the top of this pyarts-controlfile:
```
import os
print(os.environ.get("ARTS_DATA_PATH"))
```
"""
import os
import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts
# Initialize ARTS workspace
ws = pyarts.workspace.Workspace()
"""
The workspace contains an absorption species helper list
that allows it to understand what data and methods you
will want to call.
This example sets the absorption species to O2 collision-induced
absorption (CIA) with self and also with N2
"""
ws.absorption_speciesSet(species=["O2-CIA-O2", "O2-CIA-N2"])
"""
We now need to load the data.
ARTS comes with a lot of CIA data available
in its arts-cat-data repository. The first line below
downloads that data to a local cache and allows you
to run this script.
You can also provide CIA data in other ways, e.g., by loading HITRAN
or some other format. What is important is to populate the absorption
CIA data object with appropriate data.
"""
pyarts.data.download()
ws.absorption_cia_dataReadSpeciesSplitCatalog(basename="cia/")
"""
Compute absorption
Now we can use any of the methods available to compute CIA absorption.
Below we simply set a simple atmosphere and plot the resulting absorption.
Please see other examples for more details on how to use CIA
calculations in increasingly complex ways to solve your problem
"""
# Initialize an atmospheric point
atm = pyarts.arts.AtmPoint()
atm.temperature = 295 # At room temperature
atm.pressure = 1e5 # At 1 bar
atm["O2"] = 0.21 # At 21% atmospheric Oxygen
atm["N2"] = 0.78 # At 78% atmospheric Nitrogen
# Plot the absorption of this example
fig, ax = pyarts.plot(ws.absorption_cia_data, atm=atm)
for i, a in enumerate(ax.flatten()):
a.set_xlabel("Wavelength [nm]")
a.set_ylabel("Absorption [1/m]")
f0 = np.inf
f1 = -np.inf
for x in ws.absorption_cia_data[i].data:
f0 = min(f0, x.grids[0][0])
f1 = max(f1, x.grids[0][-1])
f = np.linspace(f0, f1, 7)
a.set_xticks(f, (pyarts.arts.convert.freq2wavelen(f)*1e9).round(1))
a.set_yscale('log')
a.set_title(a.get_lines()[0].get_label())
fig.suptitle("CIA Absorption")
if "ARTS_HEADLESS" not in os.environ:
plt.tight_layout()
plt.show()
(Source code, svg, pdf)
Single species recipe
"""Absorption by a single species"""
import os
import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts
# Download catalogs
pyarts.data.download()
# %% Select absorption species
species = "O2-66" # Main isotope of O2
# %% Activate the recipe for this species
#
# See [SingleSpeciesAbsorption](pyarts3.recipe.rst#pyarts3.recipe.SingleSpeciesAbsorption).
absorption = pyarts.recipe.SingleSpeciesAbsorption(species=species)
# %% Select a single temperature, a VMR value, and a range of pressures
atm = pyarts.arts.AtmPoint()
atm["O2"] = 0.2095
atm.temperature = 273
ps = np.logspace(5, -2, 8)
# %% Select frequency range
line_f0 = 118750348044.712 # Lowest energy absorption line
f = np.linspace(-500e6, 500e6, 1001) + line_f0 # Some range around it
# %% Use the recipe and convert the results to cross-sections
xsec = []
for p in ps:
atm.pressure = p
xsec.append(absorption(f, atm) / atm.number_density(species))
xsec = np.array(xsec)
# %% Plot the results
fig, ax = plt.subplots()
ax.semilogy((f - line_f0) / 1e6, xsec.T)
ax.set_xlabel("Frequency offset [MHz]")
ax.set_ylabel("Cross-section [$m^2$]")
ax.set_title("Cross-section of O$_2$ 16-16")
ax.set_ylim(ymin=1e-3 * np.min(xsec))
ax.legend(
[f"P: $10^{'{'}{round(np.log10(x))}{'}'}$" for x in ps],
ncols=4,
loc="lower center",
)
if "ARTS_HEADLESS" not in os.environ:
plt.show()
# %% Integration test by ensuring some statistics look good
assert np.isclose(6.792977548868407e-28 / xsec.mean(), 1)
assert np.isclose(5.43981642113382e-24 / xsec.sum(), 1)
assert np.isclose(1.3359834491781882e-24 / xsec.max(), 1)
assert np.isclose(2.537911691540087e-26 / xsec.std(), 1)
assert np.isclose(8.236637542411964e-35 / xsec.min(), 1)
(Source code, svg, pdf)
