Note
Go to the end to download the full example code.
Computing a simple NLTE 8542 line profile in a FAL C atmosphere
First, we import everything we need. Lightweaver is typically imported as lw, but things like the library of model atoms and Fal atmospheres need to be imported separately.
import matplotlib.pyplot as plt
import numpy as np
import lightweaver as lw
from lightweaver.fal import Falc82
from lightweaver.rh_atoms import (
Al_atom,
C_atom,
CaII_atom,
Fe_atom,
H_6_atom,
He_9_atom,
MgII_atom,
N_atom,
Na_atom,
O_atom,
S_atom,
Si_atom,
)
Now, we define the functions that will be used in our spectral synthesise. First synth_8542 which synthesises and returns the line given by an atmosphere.
def synth_8542(atmos, conserve, use_ne, wave):
"""
Synthesise a spectral line for given atmosphere with different
conditions.
Parameters
----------
atmos : lw.Atmosphere
The atmospheric model in which to synthesise the line.
conserve : bool
Whether to start from LTE electron density and conserve charge, or
simply use from the electron density present in the atomic model.
use_ne : bool
Whether to use the electron density present in the model as the
starting solution, or compute the LTE electron density.
wave : np.ndarray
Array of wavelengths over which to resynthesise the final line
profile for muz=1.
Returns
-------
ctx : lw.Context
The Context object that was used to compute the equilibrium
populations.
I_wave : np.ndarray
The intensity at muz=1 for each wavelength in `wave`.
"""
# Configure the atmospheric angular quadrature
atmos.quadrature(5)
# Configure the set of atomic models to use.
rad_set = lw.RadiativeSet(
[
H_6_atom(),
C_atom(),
O_atom(),
Si_atom(),
Al_atom(),
CaII_atom(),
Fe_atom(),
He_9_atom(),
MgII_atom(),
N_atom(),
Na_atom(),
S_atom(),
]
)
# Set H and Ca to "active" i.e. NLTE, everything else participates as an
# LTE background.
rad_set.set_active('H', 'Ca')
# Compute the necessary wavelength dependent information (SpectrumConfiguration).
spect = rad_set.compute_wavelength_grid()
# Either compute the equilibrium populations at the fixed electron density
# provided in the model, or iterate an LTE electron density and compute the
# corresponding equilibrium populations (SpeciesStateTable).
if use_ne:
eq_pops = rad_set.compute_eq_pops(atmos)
else:
eq_pops = rad_set.iterate_lte_ne_eq_pops(atmos)
# Configure the Context which holds the state of the simulation for the
# backend, and provides the python interface to the backend.
# Feel free to increase Nthreads to increase the number of threads the
# program will use.
ctx = lw.Context(atmos, spect, eq_pops, conserve_charge=conserve, Nthreads=1)
# Iterate the Context to convergence (using the iteration function now
# provided by Lightweaver)
lw.iterate_ctx_se(ctx)
# Update the background populations based on the converged solution and
# compute the final intensity for mu=1 on the provided wavelength grid.
eq_pops.update_lte_atoms_hmin_pops(atmos)
I_wave = ctx.compute_rays(wave, [atmos.muz[-1]], stokes=False)
return ctx, I_wave
The wavelength grid to output the final synthesised line on.
wave = np.linspace(853.9444, 854.9444, 1001)
Load an lw.Atmosphere object containing the FAL C atmosphere with 82 points in depth, before synthesising the Ca II 8542 AA line profile using:
The given electron density.
The electron density charge conserved from a starting LTE solution.
The LTE electron density.
These results are then plotted.
atmos_ref = Falc82()
ctx_ref, I_wave_ref = synth_8542(atmos_ref, conserve=False, use_ne=True, wave=wave)
atmos_cons = Falc82()
ctx_cons, I_wave_cons = synth_8542(atmos_cons, conserve=True, use_ne=False, wave=wave)
atmos_lte = Falc82()
ctx, I_wave_lte = synth_8542(atmos_lte, conserve=False, use_ne=False, wave=wave)
plt.plot(wave, I_wave_ref, label='Reference FAL')
plt.plot(wave, I_wave_cons, label='Reference Cons')
plt.plot(wave, I_wave_lte, label='Reference LTE n_e')
plt.legend()
plt.show()
-- Iteration 0:
dJ = 1.00e+00
(Lambda iterating background)
-- Iteration 11:
dJ = 5.41e-01
H delta = 2.3811e-01
Ca delta = 1.7571e-01
-- Iteration 22:
dJ = 1.07e-01
H delta = 5.5150e-02
Ca delta = 4.8401e-02
-- Iteration 33:
dJ = 7.59e-02
H delta = 4.6047e-02
Ca delta = 2.5978e-02
-- Iteration 44:
dJ = 6.43e-02
H delta = 3.5603e-02
Ca delta = 1.6602e-02
-- Iteration 55:
dJ = 4.86e-02
H delta = 2.4705e-02
Ca delta = 9.1835e-03
-- Iteration 66:
dJ = 2.74e-02
H delta = 1.5182e-02
Ca delta = 4.0498e-03
-- Iteration 77:
dJ = 1.18e-02
H delta = 6.7548e-03
Ca delta = 1.5101e-03
-- Iteration 88:
dJ = 4.31e-03
H delta = 2.5071e-03
Ca delta = 5.3038e-04
--------------------------------------------------------------------------------
Final Iteration: 98
--------------------------------------------------------------------------------
dJ = 1.61e-03
H delta = 9.3706e-04
Ca delta = 1.9178e-04
--------------------------------------------------------------------------------
Context converged to statistical equilibrium in 98 iterations after 1.93 s.
--------------------------------------------------------------------------------
LTE Iterations 1 (-- slowest convergence)
-- Iteration 0:
dJ = 1.00e+00
(Lambda iterating background)
-- Iteration 10:
dJ = 5.02e-01
H delta = 3.3673e-01
Ca delta = 1.7078e-01
ne delta = 1.0679e-01
-- Iteration 20:
dJ = 1.24e-01
H delta = 7.6386e-02
Ca delta = 5.3759e-02
ne delta = 1.6682e-02
-- Iteration 30:
dJ = 8.17e-02
H delta = 5.0039e-02
Ca delta = 2.6950e-02
ne delta = 1.5230e-02
-- Iteration 40:
dJ = 6.85e-02
H delta = 4.1500e-02
Ca delta = 1.5740e-02
ne delta = 1.5508e-02
-- Iteration 50:
dJ = 5.87e-02
H delta = 3.6133e-02
Ca delta = 8.3598e-03
ne delta = 1.4504e-02
-- Iteration 60:
dJ = 4.55e-02
H delta = 2.9086e-02
Ca delta = 3.6814e-03
ne delta = 1.1812e-02
-- Iteration 70:
dJ = 2.97e-02
H delta = 1.9655e-02
Ca delta = 2.2404e-03
ne delta = 7.9012e-03
-- Iteration 80:
dJ = 1.67e-02
H delta = 1.1227e-02
Ca delta = 1.4507e-03
ne delta = 4.4466e-03
-- Iteration 90:
dJ = 8.19e-03
H delta = 5.5476e-03
Ca delta = 7.5152e-04
ne delta = 2.1766e-03
-- Iteration 100:
dJ = 3.70e-03
H delta = 2.5138e-03
Ca delta = 3.5088e-04
ne delta = 9.8237e-04
-- Iteration 110:
dJ = 1.62e-03
H delta = 1.1013e-03
Ca delta = 1.5537e-04
ne delta = 4.2873e-04
--------------------------------------------------------------------------------
Final Iteration: 112
--------------------------------------------------------------------------------
dJ = 1.37e-03
H delta = 9.3152e-04
Ca delta = 1.3158e-04
ne delta = 3.6245e-04
--------------------------------------------------------------------------------
Context converged to statistical equilibrium in 112 iterations after 2.34 s.
--------------------------------------------------------------------------------
LTE Iterations 1 (-- slowest convergence)
-- Iteration 0:
dJ = 1.00e+00
(Lambda iterating background)
-- Iteration 11:
dJ = 3.52e-01
H delta = 1.8624e-01
Ca delta = 1.3708e-01
-- Iteration 22:
dJ = 1.03e-01
H delta = 5.7038e-02
Ca delta = 3.4188e-02
-- Iteration 33:
dJ = 7.19e-02
H delta = 4.7033e-02
Ca delta = 1.9780e-02
-- Iteration 44:
dJ = 5.78e-02
H delta = 3.3873e-02
Ca delta = 1.1671e-02
-- Iteration 55:
dJ = 3.90e-02
H delta = 2.0211e-02
Ca delta = 5.6154e-03
-- Iteration 66:
dJ = 1.86e-02
H delta = 1.0214e-02
Ca delta = 2.2083e-03
-- Iteration 77:
dJ = 7.23e-03
H delta = 4.0513e-03
Ca delta = 7.7384e-04
-- Iteration 86:
dJ = 3.04e-03
H delta = 1.7113e-03
Ca delta = 3.1560e-04
--------------------------------------------------------------------------------
Final Iteration: 92
--------------------------------------------------------------------------------
dJ = 1.67e-03
H delta = 9.4479e-04
Ca delta = 1.7206e-04
--------------------------------------------------------------------------------
Context converged to statistical equilibrium in 92 iterations after 1.86 s.
--------------------------------------------------------------------------------
LTE Iterations 1 (He slowest convergence)
Total running time of the script: (0 minutes 6.769 seconds)