This is perhaps less of an issue and more of a feature/docs request. My ignorance here is likely due to my unfamiliarity with Lightweaver, RH, and radiative transfer so I will apologize in advance 😅
I'm working through this toy example which uses RH (from the command line) + a combination of scripts to demonstrate the effect of the velocity on the shape of the absorption line profile around 656.25 nm. For simplicity, I'll focus just on the first part which leaves the velocity unmodified
Starting from the simple gallery example in the docs, my understanding is that I need to modify the H_6 atom that is passed to RadiativeSet in the following way:
h6 = H_6_atom()
h6.lines[4].quadrature.qCore = 25.0
h6.lines[4].quadrature.Nlambda = 200
(the 4 comes from the fact that I want to modify the first line in the Balmer series, j=2, i=1). However, my resulting line profile is
compared to the example produced in RH,
I've included the complete code to reproduce my figure below. I'm sure I'm likely doing something wrong or not completely understanding what I'm doing here!
from lightweaver.fal import Falc82
from lightweaver.rh_atoms import (H_6_atom, CaII_atom, O_atom, C_atom, Si_atom, Al_atom, Fe_atom,
He_atom, Mg_atom, N_atom, Na_atom, S_atom)
import lightweaver as lw
import matplotlib.pyplot as plt
import numpy as np
import time
def synth_8542(atmos, conserve, useNe, 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.
useNe : 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.
Iwave : 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.
# Set some things specific to our example.
h6 = H_6_atom()
h6.lines[4].quadrature.qCore = 25.0
h6.lines[4].quadrature.Nlambda = 200
aSet = lw.RadiativeSet([
h6,
C_atom(),
O_atom(),
Si_atom(),
Al_atom(),
CaII_atom(),
Fe_atom(),
He_atom(),
Mg_atom(),
N_atom(),
Na_atom(),
S_atom(),
])
# Set H and Ca to "active" i.e. NLTE, everything else participates as an
# LTE background.
aSet.set_active('H', 'Ca')
# Compute the necessary wavelength dependent information (SpectrumConfiguration).
spect = aSet.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 useNe:
eqPops = aSet.compute_eq_pops(atmos)
else:
eqPops = aSet.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, eqPops, conserveCharge=conserve, Nthreads=1)
start = time.time()
# Iterate the Context to convergence
iterate_ctx(ctx)
end = time.time()
print('%.2f s' % (end - start))
# Update the background populations based on the converged solution and
# compute the final intensity for mu=1 on the provided wavelength grid.
eqPops.update_lte_atoms_Hmin_pops(atmos)
Iwave = ctx.compute_rays(wave, [atmos.muz[-1]], stokes=False)
return ctx, Iwave
def iterate_ctx(ctx, Nscatter=3, NmaxIter=500):
'''
Iterate a Context to convergence.
'''
for i in range(NmaxIter):
# Compute the formal solution
dJ = ctx.formal_sol_gamma_matrices()
# Just update J for Nscatter iterations
if i < Nscatter:
continue
# Update the active populations under statistical equilibrium,
# conserving charge if this option was set on the Context.
delta = ctx.stat_equil()
# If we are converged in both relative change of J and populations,
# then print a message and return
# N.B. as this is just a simple case, there is no checking for failure
# to converge within the NmaxIter. This could be achieved simpy with an
# else block after this for.
if dJ < 3e-3 and delta < 1e-3:
print('%d iterations' % i)
print('-'*80)
return
atmos_rest = Falc82()
wave = np.linspace(656.1, 656.5, 1000)
ctx, Iwave_rest = synth_8542(atmos_rest, conserve=False, useNe=False, wave=wave)
plt.plot(wave, Iwave_rest)
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent