import numpy as np
from typing import Union, Tuple
try:
import QENSmodels
except ImportError:
print('Module QENSmodels not found')
[docs]def hwhmEquivalentSitesCircle(
q: Union[float, list, np.ndarray],
Nsites: int = 3,
radius: float = 1.0,
resTime: float = 1.0
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Returns some characteristics of `EquivalentSitesCircle` as functions
of the momentum transfer `q`:
the half-width half-maximum (`hwhm`), the elastic incoherent structure
factor (`eisf`), and the quasi-elastic incoherent structure factor (`qisf`)
Parameters
----------
q: float, list or :class:`~numpy:numpy.ndarray`
momentum transfer (non-fitting, in 1/Angstrom)
Nsites: integer
number of sites in circle (non-fitting). Default to 3.
radius: float
radius of the circle (in Angstrom). Default to 1.
resTime: float
residence time (in ps). Default to 1.
Returns
-------
hwhm: :class:`~numpy:numpy.ndarray`
half-width half maximum
eisf: :class:`~numpy:numpy.ndarray`
elastic incoherent structure factor
qisf: :class:`~numpy:numpy.ndarray`
quasi-elastic incoherent structure factor
Examples
--------
>>> hwhm, eisf, qisf = hwhmEquivalentSitesCircle([1., 2.], 6, 0.5, 1.5)
>>> round(hwhm[0, 1], 3)
0.333
>>> round(hwhm[0, 3], 3)
1.333
>>> round(eisf[0], 3)
0.92
>>> round(eisf[1], 3)
0.713
>>> round(qisf[0, 1], 6)
0.000503
>>> round(qisf[1, 4],6)
0.13616
"""
# input validation
q = np.asarray(q, dtype=np.float32)
if radius <= 0:
raise ValueError("radius, the radius of the circle, "
"should be positive")
if resTime < 0:
raise ValueError("resTime, the residence time, should be positive")
if Nsites < 2:
raise ValueError("the minimum number of sites N is 2")
# number of sites has to be an integer
Nsites = int(Nsites)
# index of sites in circle
sites = np.arange(Nsites)
hwhm = 2.0 / resTime * np.sin(sites * np.pi / Nsites) ** 2
hwhm = np.tile(hwhm, (q.size, 1))
# jump distances between sites
jump_distance = 2.0 * radius * np.sin(sites * np.pi / Nsites)
# QR matrix [q.size, N] and corresponding spherical Bessel functions
QR = np.outer(q, jump_distance)
sphBessel = np.ones(QR.shape)
idx = np.nonzero(QR)
sphBessel[idx] = np.sin(QR[idx]) / QR[idx]
isf = np.zeros(QR.shape)
for i in range(Nsites):
for j in range(Nsites):
isf[:, i] += sphBessel[:, j] * np.cos(2. * i * j * np.pi / Nsites)
isf[:, i] /= Nsites
eisf = isf[:, 0]
qisf = isf[:, 1:]
return hwhm, eisf, qisf
[docs]def sqwEquivalentSitesCircle(
w: Union[float, list, np.ndarray],
q: Union[float, list, np.ndarray],
scale: float = 1.0,
center: float = 0.0,
Nsites: int = 3,
radius: float = 1.0,
resTime: float = 1.0
) -> Union[float, list, np.ndarray]:
r"""
Model
`Jumps between Nsites equivalent sites on a circle with
a radius equal to `radius` `
= A_0 delta + Sum of Lorentzians ...
It models a circular random walk among these `Nsites`
sites.
Parameters
----------
w: list or :class:`~numpy:numpy.ndarray`
energy transfer (in 1/ps)
q: float, list or :class:`~numpy:numpy.ndarray`
momentum transfer (non-fitting, in 1/Angstrom)
scale: float
scale factor. Default to 1.
center: float
center of peak. Default to 0.
Nsites: integer
number of sites in circle (non-fitting). Default to 3.
radius: float
radius of rotation (in Angstrom). Default to 1.
resTime: float
residence time in a site before jumping to another site (in ps).
Default to 1.
Return
------
:class:`~numpy:numpy.ndarray`
output array
Examples
--------
>>> sqw = sqwEquivalentSitesCircle(1, 1, 1, 0, 4, 1, 1)
>>> round(sqw[0], 3)
0.045
>>> sqw = sqwEquivalentSitesCircle([1, 2, 3], [0.3, 0.4], 1, 0, 5, 1, 1)
>>> round(sqw[0, 0], 3)
0.004
>>> round(sqw[0, 1], 3)
0.001
>>> round(sqw[0, 2], 3)
0.001
>>> round(sqw[1, 0], 3)
0.008
>>> round(sqw[1, 1], 3)
0.003
>>> round(sqw[1, 2], 3)
0.001
Notes
-----
* The `sqwEquivalentSitesCircle` is expressed as
.. math::
S(q, \omega) = \text{delta}(\omega, A_0(q), \text{center} )
+ \sum_{i=1}^{N-1}
\text{Lorentzian}(\omega, A_i(q)\text{scale}, \text{center},
\Gamma_i)
where
.. math::
A_i(q) = \frac{1}{N}\sum_{j=1}^N j_0(qr_j)\cos(2ij\pi/N)
r_j = 2R \sin(j\pi/N)
\Gamma_i = \frac{2}{\text{resTime}}\sin^2(i\pi/N)
* The number of sites `N` is converted to an integer by the function.
It should **not** be used as a fitting parameter.
References
----------
* R. Hempelmann, Quasielastic Neutron Scattering and Solid State Diffusion
(Oxford, 2000).
* J. D. Barnes, *Journal of Chemical Physics* **58**, 5193-5201 (1973).
`link <https://aip.scitation.org/doi/abs/10.1063/1.1679130?journalCode=jcp>`__
""" # noqa
# Input validation
w = np.asarray(w)
q = np.asarray(q, dtype=np.float32)
# Create output array
sqw = np.zeros((q.size, w.size))
# Get widths, EISFs and QISFs of model
hwhm, eisf, qisf = hwhmEquivalentSitesCircle(q, Nsites, radius, resTime)
# Number of Lorentzians (= N-1)
numberLorentz = hwhm.shape[1] - 1
# Sum of Lorentzians
# (Note that hwhm has dimensions [q.size, N], as hwhm[:,0]
# contains a width=0, corresponding to the elastic line
# (eisf), while qisf has dimensions [q.size, N-1])
for i in range(q.size):
sqw[i, :] = eisf[i] * QENSmodels.delta(w, scale, center)
for j in range(numberLorentz):
sqw[i, :] += qisf[i, j] * QENSmodels.lorentzian(w,
scale,
center,
hwhm[i, j + 1])
# For Bumps use (needed for final plotting)
# Using a 'Curve' in bumps for each Q --> needs vector array
if q.size == 1:
sqw = np.reshape(sqw, w.size)
return sqw