"""
Decorator and function registry providing an xraylib-compatible interface backed by DABAX data.
"""
#
# dabax functions with the same interface as xraylib
#
import numpy
import scipy.constants as codata
from silx.io.specfile import SpecFile
from dabax.common_tools import atomic_symbols, atomic_number
from dabax.common_tools import bragg_metrictensor
from dabax.common_tools import calculate_f0_from_f0coeff
from dabax.dabax_base import Functions as BaseFunctions, DabaxBase
[docs]class Functions(BaseFunctions):
"""Extended function-type registry adding NIST compound data support."""
CompoundDataNIST = 'CompoundDataNIST'
[docs] @staticmethod
def get_mode(function):
"""Return the parsing mode for the given function type."""
if function == Functions.CompoundDataNIST: return 0
else: return BaseFunctions.get_mode(function)
[docs]class DabaxXraylibDecorator(object):
"""
Mixin providing an xraylib-compatible API backed by DABAX data files.
Not all xraylib functions are implemented. Mix in with :class:`DabaxBase`
(see :class:`DabaxXraylib`) to obtain a fully usable object.
"""
#########################
# crystals
#########################
[docs] def Crystal_GetCrystal(self, entry_name='YB66'):
"""
Parse a crystal structure from the DABAX file (equivalent to ``xraylib.Crystal_GetCrystal``).
Parameters
----------
entry_name : str, optional
Crystal name as stored in the DABAX Crystals file. Default ``'YB66'``.
Returns
-------
dict
Crystal dictionary with keys: ``name``, ``a``, ``b``, ``c``,
``alpha``, ``beta``, ``gamma``, ``volume``, ``n_atom``, ``atom``
(list of atom dicts each with ``Zatom``, ``fraction``, ``x``,
``y``, ``z``, ``charge``), ``cpointer``.
"""
filename = self.get_file_Crystals()
spec_file, entries = self._get_registered_spec_file(filename, Functions.Crystals)
data, labels, index_found = DabaxBase._get_data_and_labels(spec_file,
entries,
entry_name,
verbose=self.verbose(),
filename=filename)
cryst = {'name': entry_name} # returned dictionary like that one created by xraylib.Crystal_GetCrystal(descriptor)
cell_parameters = spec_file[index_found].scan_header_dict["UCELL"]
cell_parameters = ' '.join(cell_parameters.split()) # remove multiple blanks
a = cell_parameters.split(' ')
cryst['a'] = float(a[0])
cryst['b'] = float(a[1])
cryst['c'] = float(a[2])
cryst['alpha'] = float(a[3])
cryst['beta'] = float(a[4])
cryst['gamma'] = float(a[5])
volume = bragg_metrictensor(float(a[0]), float(a[1]), float(a[2]),
float(a[3]), float(a[4]), float(a[5]), RETURN_VOLUME=1)
cryst['volume'] = volume
cell_data = numpy.array(data)
cryst['n_atom'] = cell_data.shape[1]
atom = []
for i in range(cell_data.shape[1]):
if cell_data.shape[0] == 5: # standard 5 columns
# not here, this info is not in the dabax file
# s = symbol_to_from_atomic_number(int(cell_data[0,i]))
atom.append({
'Zatom': int(cell_data[0, i]),
'fraction': cell_data[1, i],
'x': cell_data[2, i],
'y': cell_data[3, i],
'z': cell_data[4, i],
'charge': 0.0, })
else: # 6 columns (charge)
# 'AtomicName' required to compatible my current code
# s = symbol_to_from_atomic_number(int(cell_data[0,i]))
# if cell_data[5, i] != 0: #charged
# s = s + f'%+.6g'%cell_data[5, i]
atom.append({
# 'AtomicName': s,
'Zatom': int(cell_data[0, i]),
'fraction': cell_data[1, i],
'x': cell_data[2, i],
'y': cell_data[3, i],
'z': cell_data[4, i],
'charge': cell_data[5, i], })
cryst['atom'] = atom
cryst['cpointer'] = None
ANISO_KEY = "UANISO_COFF" # prefix for a line with anisotropic coefficients
d = spec_file[index_found].scan_header_dict
AnisoItem = {'Name': ' ',
'start': 0,
'end': 0,
'beta11': 0.0,
'beta22': 0.0,
'beta33': 0.0,
'beta12': 0.0,
'beta13': 0.0,
'beta23': 0.0}
a = [(x, d[x].split()) for x in d if x[:len(ANISO_KEY)] == ANISO_KEY]
if len(a) > 0: # found Anisotropic coefficients in the header, process it
a = sorted(a, key=lambda x: int(x[1][0]),
reverse=False) # sort 'Start' ascendant, avoid order changed by the SpecFile
n = 0
Aniso = []
for x in a: # tuple('UANISO_COFF_B1',[1 96 0.00038 0.00044 0.00039 0 0 0])
AnisoItem['Name'] = x[0][len(ANISO_KEY) + 1:] # get site atom name starting from 13th character 'B1', etc
AnisoItem['start'] = int(x[1][0])
AnisoItem['end'] = int(x[1][1])
AnisoItem['beta11'] = float(x[1][2])
AnisoItem['beta22'] = float(x[1][3])
AnisoItem['beta33'] = float(x[1][4])
AnisoItem['beta12'] = float(x[1][5])
AnisoItem['beta13'] = float(x[1][6])
AnisoItem['beta23'] = float(x[1][7])
Aniso.append(AnisoItem.copy())
n = n + 1
cryst['Aniso'] = Aniso # if having key 'Ansio' when there is anisotropic data,otherwise no
cryst['n_aniso'] = n
else: # create a dummy Aniso to compatible with xraylib
cryst['Aniso'] = [AnisoItem.copy()]
cryst['n_aniso'] = 1
return cryst
[docs] def Crystal_GetCrystalsList(self):
"""
Return the list of crystal names available in the DABAX Crystals file.
Returns
-------
list of str
Crystal names (e.g. ``['Si', 'Ge', 'Diamond', ...]``).
"""
filename = self.get_file_Crystals()
_, crystals = self._get_registered_spec_file(filename, Functions.Crystals)
return crystals
[docs] def Crystal_dSpacing(self, cryst, h, k, l):
"""
Return the d-spacing for Miller indices (h, k, l).
Parameters
----------
cryst : dict
Crystal dictionary as returned by :meth:`Crystal_GetCrystal`.
h, k, l : int
Miller indices.
Returns
-------
float
d-spacing in Å.
"""
return bragg_metrictensor(cryst['a'], cryst['b'], cryst['c'],
cryst['alpha'], cryst['beta'], cryst['gamma'],
HKL=[h, k, l])
[docs] def Bragg_angle(self, cryst, E_keV, h, k, l):
"""
Return the Bragg angle for a given crystal reflection and photon energy.
Parameters
----------
cryst : dict
Crystal dictionary as returned by :meth:`Crystal_GetCrystal`.
E_keV : float
Photon energy in keV.
h, k, l : int
Miller indices.
Returns
-------
float
Bragg angle in radians.
"""
dspacing = self.Crystal_dSpacing(cryst, h, k, l) # in A
wavelength = codata.h * codata.c / codata.e / (E_keV * 1e3) * 1e10 # in A
return numpy.arcsin(wavelength / 2 / dspacing)
[docs] def Crystal_F_H_StructureFactor(self,
crystal_id,
energy_in_kev,
millerH,
millerK,
millerL,
debyeWaller,
ratio_theta_thetaB=1.0):
"""
Calculate the crystal structure factor F_H for a given reflection.
Parameters
----------
crystal_id : dict
Crystal dictionary as returned by :meth:`Crystal_GetCrystal`.
energy_in_kev : float
Photon energy in keV.
millerH, millerK, millerL : int
Miller indices.
debyeWaller : float
Debye-Waller factor.
ratio_theta_thetaB : float, optional
Ratio of theta to the Bragg angle. Default 1.0.
Returns
-------
complex or numpy.ndarray
Structure factor F_H.
"""
energy = energy_in_kev * 1e3
wavelength = codata.h * codata.c / codata.e / energy * 1e10
# print(crystal_id["n_atom"])
atom = crystal_id['atom']
natom = len(atom)
list_fraction = [atom[i]['fraction'] for i in range(natom)]
list_x = [atom[i]['x'] for i in range(natom)]
list_y = [atom[i]['y'] for i in range(natom)]
list_z = [atom[i]['z'] for i in range(natom)]
F_H = numpy.zeros(numpy.array(energy).size, dtype=complex)
for i in range(natom):
atom_i = atom[i]
if (i > 0) and (atom_i['Zatom'] == Z_i) and (atom_i['charge'] == charge_i):
pass # avoid re-calculating f0 if inputs are identical to previous call
else:
Z_i = atom_i['Zatom']
charge_i = atom_i['charge']
coeffs = self.f0_with_fractional_charge(Z_i, charge=charge_i)
if (millerH == 0 and millerK == 0 and millerL == 0):
ratio = 0.0
else:
angle = self.Bragg_angle(crystal_id, energy_in_kev,
millerH, millerK, millerL) * ratio_theta_thetaB
ratio = numpy.sin(angle) / wavelength
f0_i = calculate_f0_from_f0coeff(coeffs, ratio)
Fi = self.Fi(Z_i, energy_in_kev)
Fii = self.Fii(Z_i, energy_in_kev)
F_H += (f0_i + Fi - Fii * 1j) * list_fraction[i] * \
numpy.exp(2j*numpy.pi*(millerH*list_x[i]+millerK*list_y[i]+millerL*list_z[i]))
return F_H * debyeWaller
#########################
# misc
#########################
[docs] def CompoundParser(self, descriptor):
"""
Parse a chemical formula string into its elemental composition.
Parameters
----------
descriptor : str
Chemical formula string (e.g. ``"H2O"``, ``"Ca5(PO4)3F"``).
Returns
-------
dict
Dictionary with keys ``nElements``, ``nAtomsAll``, ``Elements``
(list of atomic numbers Z), ``massFractions``, ``nAtoms``,
``molarMass``.
"""
return self.compound_parser(descriptor)
[docs] def SymbolToAtomicNumber(self, symbol):
"""
Return the atomic number corresponding to an element symbol.
Parameters
----------
symbol : str
Element symbol (e.g. ``"Si"``).
Returns
-------
int
Atomic number Z.
"""
return atomic_number(symbol)
[docs] def AtomicNumberToSymbol(self, Z):
"""
Return the element symbol for a given atomic number.
Parameters
----------
Z : int
Atomic number.
Returns
-------
str
Element symbol (e.g. ``"Si"`` for Z=14).
"""
return atomic_symbols()[Z]
[docs] def ElementDensity(self, Z):
"""
Return the density of an element at room temperature.
Parameters
----------
Z : int
Atomic number.
Returns
-------
float
Density in g/cm³.
"""
return self.element_density(self.AtomicNumberToSymbol(Z))
[docs] def AtomicWeight(self, Z):
"""
Return the atomic weight of an element.
Parameters
----------
Z : int
Atomic number.
Returns
-------
float
Atomic weight in g/mol.
"""
return self.atomic_weights(self.AtomicNumberToSymbol(Z))
#########################
# scattering functions
#########################
[docs] def Fi(self, Z, energy):
"""
Return the real part of the anomalous scattering factor correction Δf′.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Δf′ (real anomalous dispersion correction).
"""
return self.FiAndFii(Z, energy)[0]
[docs] def Fii(self, Z, energy):
"""
Return the imaginary part of the anomalous scattering factor correction Δf″.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Δf″ (imaginary anomalous absorption correction, always positive).
"""
return self.FiAndFii(Z, energy)[1]
[docs] def FF_Rayl(self, Z, q):
"""
Return the atomic form factor for Rayleigh (coherent) scattering.
Parameters
----------
Z : int
Atomic number.
q : float or numpy.ndarray
Momentum transfer in Å⁻¹ (q = sin(θ)/λ).
Returns
-------
float or numpy.ndarray
Atomic form factor f0(q).
"""
coeffs = self.f0_with_fractional_charge(Z, charge=0.0)
return calculate_f0_from_f0coeff(coeffs, q)
#########################
# cross sections
#########################
# barn/atom
[docs] def CSb_Total(self, Z, energy):
"""
Return the total attenuation cross section in barn/atom.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Total cross section in barn/atom.
"""
return self.crosssec_interpolate(self.AtomicNumberToSymbol(Z), energy * 1e3,
partial='TotalCrossSection[barn/atom]',)
[docs] def CSb_Photo(self, Z, energy):
"""
Return the photoionization cross section in barn/atom.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Photoionization cross section in barn/atom.
"""
return self.crosssec_interpolate(self.AtomicNumberToSymbol(Z), energy * 1e3,
partial='PhotoElectric[barn/atom]',)
[docs] def CSb_Rayl(self, Z, energy):
"""
Return the Rayleigh (coherent) scattering cross section in barn/atom.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Rayleigh scattering cross section in barn/atom.
"""
return self.crosssec_interpolate(self.AtomicNumberToSymbol(Z), energy * 1e3,
partial='Rayleigh(coherent)[barn/atom]',)
[docs] def CSb_Compt(self, Z, energy):
"""
Return the Compton (incoherent) scattering cross section in barn/atom.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Compton scattering cross section in barn/atom.
"""
return self.crosssec_interpolate(self.AtomicNumberToSymbol(Z), energy * 1e3,
partial='Compton(incoherent)[barn/atom]',)
# cm²/g (mass attenuation)
[docs] def CS_Total(self, Z, energy):
"""
Return the total mass attenuation cross section in cm²/g.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Total mass attenuation cross section in cm²/g.
"""
return self.CSb_Total(Z, energy) * (codata.Avogadro * 1e-24 / self.AtomicWeight(Z))
[docs] def CS_Photo(self, Z, energy):
"""
Return the photoionization mass cross section in cm²/g.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Photoionization mass cross section in cm²/g.
"""
return self.CSb_Photo(Z, energy) * (codata.Avogadro * 1e-24 / self.AtomicWeight(Z))
[docs] def CS_Rayl(self, Z, energy):
"""
Return the Rayleigh (coherent) scattering mass cross section in cm²/g.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Rayleigh scattering mass cross section in cm²/g.
"""
return self.CSb_Rayl(Z, energy) * (codata.Avogadro * 1e-24 / self.AtomicWeight(Z))
[docs] def CS_Compt(self, Z, energy):
"""
Return the Compton (incoherent) scattering mass cross section in cm²/g.
Parameters
----------
Z : int
Atomic number.
energy : float
Photon energy in keV.
Returns
-------
float
Compton scattering mass cross section in cm²/g.
"""
return self.CSb_Compt(Z, energy) * (codata.Avogadro * 1e-24 / self.AtomicWeight(Z))
# for compounds — descriptor can be a formula string or a NIST compound name
[docs] def CS_Total_CP(self, descriptor, energy):
"""
Return the total mass attenuation cross section in cm²/g for a compound.
Parameters
----------
descriptor : str
Chemical formula (e.g. ``"SiO2"``) or NIST compound name.
energy : float
Photon energy in keV.
Returns
-------
float
Total mass attenuation cross section in cm²/g.
"""
cp = self.CompoundParserCheckingNIST(descriptor,)
out = 0.0
for i in range(cp["nElements"]):
out += self.CS_Total(cp["Elements"][i], energy) * cp["massFractions"][i]
return out
[docs] def CSb_Total_CP(self, descriptor, energy):
"""
Return the total attenuation cross section in barn/atom for a compound.
Parameters
----------
descriptor : str
Chemical formula or NIST compound name.
energy : float
Photon energy in keV.
Returns
-------
float
Total cross section in barn/atom.
"""
cp = self.CompoundParserCheckingNIST(descriptor,)
out = 0.0
for i in range(cp["nElements"]):
out += self.CSb_Total(cp["Elements"][i], energy) * cp["massFractions"][i]
return out
[docs] def CS_Photo_CP(self, descriptor, energy):
"""
Return the photoionization mass cross section in cm²/g for a compound.
Parameters
----------
descriptor : str
Chemical formula or NIST compound name.
energy : float
Photon energy in keV.
Returns
-------
float
Photoionization mass cross section in cm²/g.
"""
cp = self.CompoundParserCheckingNIST(descriptor,)
out = 0.0
for i in range(cp["nElements"]):
out += self.CS_Photo(cp["Elements"][i], energy) * cp["massFractions"][i]
return out
[docs] def CSb_Photo_CP(self, descriptor, energy):
"""
Return the photoionization cross section in barn/atom for a compound.
Parameters
----------
descriptor : str
Chemical formula or NIST compound name.
energy : float
Photon energy in keV.
Returns
-------
float
Photoionization cross section in barn/atom.
"""
cp = self.CompoundParserCheckingNIST(descriptor,)
out = 0.0
for i in range(cp["nElements"]):
out += self.CSb_Photo(cp["Elements"][i], energy) * cp["massFractions"][i]
return out
[docs] def CS_Rayl_CP(self, descriptor, energy):
"""
Return the Rayleigh scattering mass cross section in cm²/g for a compound.
Parameters
----------
descriptor : str
Chemical formula or NIST compound name.
energy : float
Photon energy in keV.
Returns
-------
float
Rayleigh scattering mass cross section in cm²/g.
"""
cp = self.CompoundParserCheckingNIST(descriptor,)
out = 0.0
for i in range(cp["nElements"]):
out += self.CS_Rayl(cp["Elements"][i], energy) * cp["massFractions"][i]
return out
[docs] def CSb_Rayl_CP(self, descriptor, energy):
"""
Return the Rayleigh scattering cross section in barn/atom for a compound.
Parameters
----------
descriptor : str
Chemical formula or NIST compound name.
energy : float
Photon energy in keV.
Returns
-------
float
Rayleigh scattering cross section in barn/atom.
"""
cp = self.CompoundParserCheckingNIST(descriptor,)
out = 0.0
for i in range(cp["nElements"]):
out += self.CSb_Rayl(cp["Elements"][i], energy) * cp["massFractions"][i]
return out
[docs] def CS_Compt_CP(self, descriptor, energy):
"""
Return the Compton scattering mass cross section in cm²/g for a compound.
Parameters
----------
descriptor : str
Chemical formula or NIST compound name.
energy : float
Photon energy in keV.
Returns
-------
float
Compton scattering mass cross section in cm²/g.
"""
cp = self.CompoundParserCheckingNIST(descriptor,)
out = 0.0
for i in range(cp["nElements"]):
out += self.CS_Compt(cp["Elements"][i], energy) * cp["massFractions"][i]
return out
[docs] def CSb_Compt_CP(self, descriptor, energy):
"""
Return the Compton scattering cross section in barn/atom for a compound.
Parameters
----------
descriptor : str
Chemical formula or NIST compound name.
energy : float
Photon energy in keV.
Returns
-------
float
Compton scattering cross section in barn/atom.
"""
cp = self.CompoundParserCheckingNIST(descriptor,)
out = 0.0
for i in range(cp["nElements"]):
out += self.CSb_Compt(cp["Elements"][i], energy) * cp["massFractions"][i]
return out
#########################
# refractive index
#########################
[docs] def Refractive_Index(self, descriptor, energy, density):
"""
Return the complex refractive index n = n_re + i·n_im for a material.
Parameters
----------
descriptor : str
Chemical formula (e.g. ``"Be"``, ``"SiO2"``) or NIST compound name.
energy : float
Photon energy in keV.
density : float
Material density in g/cm³.
Returns
-------
complex
Complex refractive index n.
"""
cp = self.compound_parser(descriptor)
KD = 4.15179082788e-4 # TODO: recalculate with codata....
rv_re = 0.0
rv_im = 0.0
for i in range(cp["nElements"]):
Z = cp["Elements"][i]
rv_re += cp["massFractions"][i] * KD * (Z + self.Fi(Z, energy)) / self.AtomicWeight(Z) / energy / energy
rv_im += self.CS_Total(Z, energy) * cp["massFractions"][i]
return (1 - rv_re * density) + 1j*(rv_im * density * 9.8663479e-9 / energy)
[docs] def Refractive_Index_Re(self, descriptor, energy, density):
"""
Return the real part of the refractive index for a material.
Parameters
----------
descriptor : str
Chemical formula or NIST compound name.
energy : float
Photon energy in keV.
density : float
Material density in g/cm³.
Returns
-------
float
Real part of the refractive index (≤ 1 for X-rays).
"""
cp = self.compound_parser(descriptor)
KD = 4.15179082788e-4 # TODO: recalculate with codata....
rv = 0.0
for i in range(cp["nElements"]):
Z = cp["Elements"][i]
rv += cp["massFractions"][i] * KD * (Z + self.Fi(Z, energy)) / \
self.AtomicWeight(Z) / energy / energy
return (1 - rv * density)
[docs] def Refractive_Index_Im(self, descriptor, energy, density):
"""
Return the imaginary part of the refractive index for a material.
Parameters
----------
descriptor : str
Chemical formula or NIST compound name.
energy : float
Photon energy in keV.
density : float
Material density in g/cm³.
Returns
-------
float
Imaginary part of the refractive index (related to absorption).
"""
cp = self.compound_parser(descriptor)
rv = 0.0
for i in range(cp["nElements"]):
Z = cp["Elements"][i]
rv += self.CS_Total(Z, energy) * cp["massFractions"][i]
# /*9.8663479e-9 is calculated as planck's constant * speed of light / 4Pi */
# return rv * density * 9.8663479e-9 / E;
return rv * density * 9.8663479e-9 / energy
#
# NIST compounds
#
[docs] def GetCompoundDataNISTList(self):
"""
Return the list of NIST compound names available in the DABAX file.
Returns
-------
list of str
NIST compound names (e.g. ``['Air, Dry', 'Water, Liquid', ...]``).
"""
filename = self.get_file_NIST()
_, compound_list = self._get_registered_spec_file(filename, Functions.CompoundDataNIST)
return compound_list
[docs] def GetCompoundDataNISTByIndex(self, index_found):
"""
Return NIST compound composition and density data by catalog index.
Parameters
----------
index_found : int
Zero-based index into the NIST compound list.
Returns
-------
dict
Dictionary with keys: ``name`` (str), ``nElements`` (int),
``density`` (float, g/cm³), ``Elements`` (list of Z),
``massFractions`` (list of float).
"""
filename = self.get_file_NIST()
spec_file, _ = self._get_registered_spec_file(filename, Functions.CompoundDataNIST)
s1 = spec_file[index_found]
data = s1.data
name = s1.scan_header_dict["Uname"]
nElements = int(s1.scan_header_dict["UnElements"])
density = float(s1.scan_header_dict["Udensity"])
Elements = []
massFractions = []
for i in range(nElements):
Elements.append(int(data[0][i]))
massFractions.append(data[1][i])
return {"name": name, 'nElements': nElements, 'density': density, 'Elements': Elements, 'massFractions': massFractions}
[docs] def GetCompoundDataNISTByName(self, entry_name):
"""
Return NIST compound composition and density data by compound name.
Parameters
----------
entry_name : str
NIST compound name (e.g. ``'Water, Liquid'``).
Returns
-------
dict
Same structure as :meth:`GetCompoundDataNISTByIndex`.
"""
return self.GetCompoundDataNISTByIndex(self.GetCompoundDataNISTList().index(entry_name))
#
#
# DONE:
#
#
# (used in xoppy_xraylib_util):
# xraylib.Crystal_GetCrystal(descriptor)
# xraylib.Crystal_dSpacing(cryst, hh, kk, ll)
# xraylib.Crystal_dSpacing
# xraylib.CompoundParser(descriptor)
# xraylib.SymbolToAtomicNumber(descriptor)
# xraylib.AtomicNumberToSymbol(zi)
# xraylib.ElementDensity(Z)
# xraylib.AtomicWeight
# xraylib.FF_Rayl(xraylib.SymbolToAtomicNumber(descriptor), iqscale)
# xraylib.Fi(Z,1e-3*ienergy)
# xraylib.Fii(Z,1e-3*ienergy)
# xraylib.Crystal_F_H_StructureFactor(_crystal, E_keV, h, k, l, _debyeWaller, 1.0)
# xraylib.Crystal_F_H_StructureFactor(_crystal, E_keV, h, k, l, _debyeWaller, 1.0)
#
# (used in power/power3d)
#
# xraylib.CS_Total(Z,1e-3*ienergy)
# xraylib.CSb_Total(Z,1e-3*ienergy)
# xraylib.CS_Total_CP(Z,1e-3*ienergy)
# xraylib.CSb_Total_CP(Z,1e-3*ienergy)
# xraylib.Refractive_Index_Re(descriptor, energy_in_keV, density)
# xraylib.Refractive_Index_Im(descriptor, energy_in_keV, density)
#
# (used in calc_cross_sec)
# xraylib.CS_Phot()
# xraylib.CSb_Photo(Z,1e-3*ienergy)
# xraylib.CS_Phot_CP()
# xraylib.CSb_Photo_CP(descriptor,1e-3*ienergy)
# xraylib.CSb_Photo_CP(descriptor,1e-3*ienergy)
# xraylib.CS_Rayl(Z,1e-3*ienergy)
# xraylib.CSb_Rayl(Z,1e-3*ienergy)
# xraylib.CS_Rayl_CP(descriptor,1e-3*ienergy)
# xraylib.CSb_Rayl_CP(descriptor,1e-3*ienergy)
# xraylib.CS_Compt(Z,1e-3*ienergy)
# xraylib.CSb_Compt(Z,1e-3*ienergy)
# xraylib.CS_Compt_CP(descriptor,1e-3*ienergy)
# xraylib.CSb_Compt_CP(descriptor,1e-3*ienergy)
# xraylib.GetCompoundDataNISTList()
# xraylib.GetCompoundDataNISTByName(DESCRIPTOR)
# xraylib.GetCompoundDataNISTByIndex(DESCRIPTOR)
#
# TODO
#
#
# auxiliar methods
# there are not in xraylib, but accelerate the calculation
#
[docs] def CompoundParserCheckingNIST(self, descriptor):
"""
Parse a compound descriptor, trying formula parsing first, then the NIST database.
Parameters
----------
descriptor : str
Chemical formula (e.g. ``"H2O"``) or NIST compound name
(e.g. ``"Water, Liquid"``).
Returns
-------
dict
Composition dictionary as returned by :meth:`CompoundParser` or
:meth:`GetCompoundDataNISTByName`.
"""
try:
out_dict = self.compound_parser(descriptor)
except:
try:
out_dict = self.GetCompoundDataNISTByName(descriptor)
except:
raise Exception("Error processing compound descriptor: %s" % descriptor)
return out_dict
[docs] def FiAndFii(self, Z, energy):
"""
Return both anomalous scattering factor corrections (Δf′, Δf″) in one call.
Parameters
----------
Z : int
Atomic number.
energy : float or numpy.ndarray
Photon energy in keV.
Returns
-------
tuple of float or numpy.ndarray
``(f1, f2)`` where ``f1`` = Δf′ and ``f2`` = Δf″ (both positive).
"""
symbol = self.AtomicNumberToSymbol(Z)
f1, f2 = self.f1f2_interpolate(symbol, energy*1e3)
if self.get_file_f1f2() in ['f1f2_Windt.dat','f1f2_Henke.dat','f1f2_EPDL97.dat','f1f2_Chantler.dat','f1f2_asf_Kissel.dat']:
f1 -= Z
f2 *= -1.0
return f1,f2
[docs] def Crystal_F_0_F_H_F_H_bar_StructureFactor(self,
crystal_id,
energy_in_kev,
millerH,
millerK,
millerL,
debyeWaller,
rel_angle=1.0):
"""
Calculate F_0, F_H, and F_H_bar structure factors in one call.
Parameters
----------
crystal_id : dict
Crystal dictionary as returned by :meth:`Crystal_GetCrystal`.
energy_in_kev : float
Photon energy in keV.
millerH, millerK, millerL : int
Miller indices.
debyeWaller : float
Debye-Waller factor.
rel_angle : float, optional
Ratio of theta to the Bragg angle. Default 1.0.
Returns
-------
tuple of complex
(F_0, F_H, F_H_bar) structure factors.
"""
energy = energy_in_kev * 1e3
wavelength = codata.h * codata.c / codata.e / energy * 1e10
# print(crystal_id["n_atom"])
atom = crystal_id['atom']
natom = len(atom)
list_fraction = [atom[i]['fraction'] for i in range(natom)]
list_x = [atom[i]['x'] for i in range(natom)]
list_y = [atom[i]['y'] for i in range(natom)]
list_z = [atom[i]['z'] for i in range(natom)]
F_0 = numpy.zeros(numpy.array(energy).size, dtype=complex)
F_H = numpy.zeros(numpy.array(energy).size, dtype=complex)
F_H_bar = numpy.zeros(numpy.array(energy).size, dtype=complex)
for i in range(natom):
atom_i = atom[i]
if (i > 0) and (atom_i['Zatom'] == Z_i) and (atom_i['charge'] == charge_i):
pass # avoid re-calculating f0 if inputs are identical to previous call
else:
Z_i = atom_i['Zatom']
charge_i = atom_i['charge']
coeffs = self.f0_with_fractional_charge(Z_i, charge=charge_i)
if (millerH == 0 and millerK == 0 and millerL == 0):
ratio = 0.0
else:
angle = self.Bragg_angle(crystal_id, energy_in_kev,
millerH, millerK, millerL)
ratio = numpy.sin(angle * rel_angle) / wavelength
f0_i_zero = calculate_f0_from_f0coeff(coeffs, 0.0)
f0_i = calculate_f0_from_f0coeff(coeffs, ratio)
Fi = self.Fi(Z_i, energy_in_kev)
Fii = self.Fii(Z_i, energy_in_kev)
F_0 += (f0_i_zero + Fi - Fii * 1j) * list_fraction[i] * debyeWaller
F_H += (f0_i + Fi - Fii * 1j) * list_fraction[i] * debyeWaller * \
numpy.exp(2j*numpy.pi*(millerH*list_x[i]+millerK*list_y[i]+millerL*list_z[i]))
F_H_bar += (f0_i + Fi - Fii * 1j) * list_fraction[i] * debyeWaller * \
numpy.exp(2j*numpy.pi*(-millerH*list_x[i]-millerK*list_y[i]-millerL*list_z[i]))
return F_0, F_H, F_H_bar