-
Notifications
You must be signed in to change notification settings - Fork 25
Tutorial_XRTS
In this example, we demonstrate how to use the PlasmaXRTSCalculator and related classes to setup, execute, and diagnose simulation of XFEL photons scattering inelastically from a homogeneous plasma characterized by the chemical composition of the plasma, the average electron density, the average charge state, the electron temperature, and the ion temperature. The backengine is the code "xrs", a copy of which can be obtained free of charge from the maintainers of simex_platform.
In our example, we set up a calculation for a Be plasma with an average charge state 2, isochorically heated (i.e. no compression or expansion with respect to the solid state) to an electron temperature of 13 eV and ion temperature of 6 eV. XFEL photons of 4.96 keV scatter from the plasma and the scattered radiation is collected in an energy dispersive detector.
The PlasmaXRTSCalculator accepts wavefront data in native WPG format (i.e. not openPMD) as input to extract the photon energy spectrum. The total scattering spectrum will then be calculated as the convolution of the dynamic structure factor (DSF) of the plasma and the source spectrum. The former will be calculated in a semi-analytic way. A variety of approximations to the DSF are implemented in xrs, which are exposed to the user via the PlasmaXRTSCalculatorParameters. All parameters gouverning the calculation of the DSF are members of this class.
We first import the required python modules:
# SIMEX example script for calculation of X-ray Thomson Scattering from a plasma
from SimEx import PlasmaXRTSCalculatorParameters
from SimEx import PlasmaXRTSCalculator
Information about the scattering photons are read in from a wavefront file, e.g. one calculated i with the XFELPhotonPropagator and stored in native WPG hdf5 format.
source_input = "prop_out_0000001.h5"
The PlasmaXRTSCalculator depends on many parameters that characterize the physical state of the plasma, or tell the backengine which approximation to apply to calculate the DSF. Consult the reference manual section about the PlasmaXRTSCalculatorParameters for a detailed explanation of all options. In this tutorial, we will use the Born-Mermin Approximation for the high frequency (plasmon) part of the DSF, the Rayleigh peak will be calculated using the Debye-Hueckel theory and the bound-free part will be calculated from the impulse approximation.
parameters = PlasmaXRTSCalculatorParameters(
elements=[['Be', 1, -1]], # Stochiometry and partial charges
photon_energy=4960.0, # [eV]
scattering_angle=30.0, # [deg]
electron_temperature=13.0, # [eV/kB]
electron_density=3.0e23, # [1/cm**3]
ion_temperature=6.0, # [eV]
ion_charge=2.0,
mass_density=1.85, # [g/cm**3]
debye_temperature=None,
band_gap=None,
energy_range={'min' : -200.0, # Min. energy/eV to calculate (relative to photon energy)
'max' : 200.0, # Max. energy/eV to calculate (relative to photon energy)
'step': 1.0}, # Energy binning/eV.
model_Sii='DH', # Use Debye-Hueckel
model_See='BMA', # Use Born-Mermin
model_Sbf='IA', # Use impulse approximation
model_IPL=0.0, # No ionization potential lowering.
model_Mix=None, # Use default (advanced mixing).
lfc=None, # No local field correction.
Sbf_norm=None, # No normalization of the bound-free spectrum.
source_spectrum='GAUSS', # Source spectrum will be taken from wavefront input.
source_spectrum_fwhm=5, # Not needed here.
)
We construct the calculator and execute the backengine method:
xrts_calculator = PlasmaXRTSCalculator(parameters=parameters,
input_path=source_input,
output_path='Be_xrts.h5')
xrts_calculator._readH5() # Loads the source spectrum
xrts_calculator.backengine()
xrts_calculator.saveH5()
The last command saves the data in hdf5 format in the location specified by the output_path parameter of PlasmaXRTSCalculator.
Instead of re-loading the data from the hdf5 file into memory and then plotting it, we can also analyse the data directly after the backengine() command has returned:
data = xrts_calculator.data
energies = data[:,0]
See = data[:,1]
Sbf = data[:,2]
Stot = data[:,3]
import pylab
pylab.plot(energies, See, label="free-free")
pylab.plot(energies, Sbf, label="bound-free")
pylab.plot(energies, Stot, label="total")
pylab.xlabel("energy (eV)")
pylab.ylabel(r"$S(k,\omega)$ (1/eV)")
pylab.legend()
The resulting figure will look similar to this one:
Finally, for the impatient, here's the link to the entire script .