# -*- coding: utf-8 -*-
# written by Ralf Biehl at the Forschungszentrum Jülich ,
# Jülich Center for Neutron Science (JCNS-1)
# Jscatter is a program to read, analyse and plot data
# Copyright (C) 2015-2026 Ralf Biehl
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import inspect
import math
import os
import sys
import numbers
import numpy as np
from numpy import linalg as la
import scipy
import scipy.integrate
import scipy.interpolate
import scipy.constants
import scipy.special as special
import scipy.integrate as integrate
try:
from scipy.special import sph_harm_y as Ylm
except ImportError:
import scipy.special
def Ylm(n, m, theta, phi):
return scipy.special.sph_harm(m, n, theta, phi)
from .. import dataArray as dA
from .. import dataList as dL
from .. import formel
from .. formel import convolve, lorentz
from .fft import k
try:
from ..libs import fscatter
useFortran = True
except ImportError:
useFortran = False
__all__ = ['elastic_w', 'transDiff_w', 'jumpDiff_w',
'diffusionHarmonicPotential_w', 'diffusionInSphere_w',
'rotDiffusion_w', 'nSiteJumpDiffusion_w', 'resolution_w', 'lorentz_w', 'stretchedExp_w',
'doubleStretchedExp_w', 'threeLorentz_w','fixedWindowScanArrhenius','stretchedExp_beta05_w']
pi = np.pi
_path_ = os.path.realpath(os.path.dirname(__file__))
#: Planck constant in µeV*ns
h = scipy.constants.Planck / scipy.constants.e * 1E15 # µeV*ns
#: h/2π reduced Planck constant in µeV*ns
hbar = h/2/pi # µeV*ns
try:
# change in scipy 18
spjn = special.spherical_jn
except AttributeError:
spjn = lambda n, z: special.jv(n + 1 / 2, z) * np.sqrt(pi / 2) / (np.sqrt(z))
[docs]
def elastic_w(w):
r"""
Elastic line; dynamic structure factor in w domain.
:math:`\delta(w)` distribution at I(w)=0 except I(0)=a that np.trapezoid(I(w)) =1
Parameters
----------
w : array
Frequencies in 1/ns.
Zero value should be part of w.
Notes
-----
For shifted frequencies (e.g. to be symmetric around zero)
the elastic peak position (zero value) might not appear on the X values.
Returns
-------
dataArray
"""
whereis0 = np.isclose(w, 0)
if not np.any(whereis0):
raise ValueError('No zero in frequencies. elastic_w requires zero in frequencies!')
i = np.argmax(whereis0)
Iqw = np.zeros_like(w)
Iqw[i] = 1./ np.diff(w[i-1:i+2]).mean() # that we get np.trapezoid() = 1
result = dA(np.c_[w, Iqw].T)
result.setColumnIndex(iey=None)
result.columnname = 'w;Iqw'
result.modelname = inspect.currentframe().f_code.co_name
return result
[docs]
def transDiff_w(w, q, D):
r"""
Translational diffusion; dynamic structure factor in w domain.
Parameters
----------
w : array
Frequencies in 1/ns
q : float
Wavevector in nm**-1
D : float
Diffusion constant in nm**2/ns
Returns
-------
dataArray
Notes
-----
Equ 33 in [1]_
.. math:: I(\omega,q) = \frac{1}{\pi} \frac{Dq^2}{(Dq^2)^2 + \omega^2}
Examples
--------
::
import jscatter as js
import numpy as np
ql = np.r_[0.5,2:18.:3]
w = np.r_[-100:100:0.1]
p = js.grace()
iqwD = js.dL([js.dynamic.transDiff_w(w,q,0.02) for q in ql])
p.plot(iqwD,le=f'q=$wavevector nm\S-1')
p.yaxis(scale='l',label=r'S(\xw\f{}) / a.u.',min=1e-7,max=100)
p.xaxis(label=r'\xw\f{} / ns\S-1',min=-100,max=100)
p.legend(x=30,y=100,charsize=0.8)
# p.save(js.examples.imagepath+'/transDiff_w.jpg', size=(2,2))
.. image:: ../../examples/images/transDiff_w.jpg
:align: center
:width: 50 %
:alt: dynamic_t2f_examples.
References
----------
.. [1] Scattering of Slow Neutrons by a Liquid
Vineyard G Physical Review 1958 vol: 110 (5) pp: 999-1010
"""
dw = q * q * D
res = 1 / pi * dw / (dw * dw + w * w)
result = dA(np.c_[w, res].T)
result.setColumnIndex(iey=None)
result.columnname = 'w;Iqw'
result.modelname = inspect.currentframe().f_code.co_name
result.wavevector = q
result.D = D
return result
[docs]
def jumpDiff_w(w, q, t0, r0):
r"""
Jump diffusion; dynamic structure factor in w domain.
Jump diffusion as a Markovian random walk. Jump length distribution is a Gaussian
with width r0 and jump rate distribution with width G (Poisson).
Diffusion coefficient D=r0**2/2t0.
Parameters
----------
w : array
Frequencies in 1/ns
q : float
Wavevector in nm**-1
t0 : float
Mean residence time in a Poisson distribution of jump times. In units ns.
G = 1/tg = Mean jump rate
r0 : float
Root mean square jump length in 3 dimensions <r**2> = 3*r_0**2
Returns
-------
dataArray
Notes
-----
Equ 6 + 8 in [1]_ :
.. math:: S_{inc}(q,\omega) = \frac{1}{\pi} \frac{\Delta\omega}{\Delta\omega^2 + \omega^2}
with \; \Delta\omega = \frac{1-e^{-q^2 r_0^2/2}}{t_0}
References
----------
.. [1] Incoherent neutron scattering functions for random jump diffusion in bounded and infinite media.
Hall, P. L. & Ross, D. K. Mol. Phys. 42, 637–682 (1981).
"""
Ln = lambda w, dw: dw / (dw * dw + w * w) / pi
dw = 1. / t0 * (1 - np.exp(-q ** 2 * r0 ** 2 / 2.))
result = dA(np.c_[w, Ln(w, dw)].T)
result.modelname = inspect.currentframe().f_code.co_name
result.setColumnIndex(iey=None)
result.columnname = 'w;Iqw'
result.wavevector = q
result.meanresidencetime = t0
result.meanjumplength = r0
return result
[docs]
def stretchedExp_beta05_w(w, gamma):
r"""
KWW in frequency domain as analytic Fourier tranform of stretched exponential with beta=0.5.
Same as Kohlrausch--Williams--Watts (KWW) in frequency domain as
Fourier tranform of stretched exponential with beta=1.
We use here the normalized form :math:`\int_{-\infty}^{\infty} \tilde{\phi}(\omega) = 1`
Parameters
----------
w : float or array-like
Frequency
gamma : float
Relaxation rate :math:`\gamma` in units 1/[unit t] with characteristic time :math:`\tau=1/\gamma` .
Returns
-------
float or ndarray, always >= 0
Notes
-----
Analytical derivation of the **cosine Fourier transform of the Kohlrausch--Williams--Watts (KWW)**
stretched exponential with :math:`\beta = 1/2`, expressed in closed form via Fresnel integrals.
The solution can be found in [deGennes19]_ equ. 11.
We write here in a normalized form :math:`\int_{-\infty}^{\infty} \tilde{\phi}(\omega) = 1`
Below the derivation is given with missing prefactor :math:`\frac{1}{\pi}` :
.. math::
\tilde{\phi}(\omega)
= \frac{1}{\pi} \sqrt{\frac{\pi}{2\omega^3\tau}}
\left[
\left(\tfrac{1}{2} - S(\theta)\right)\sin\theta
+
\left(\tfrac{1}{2} - C(\theta)\right)\cos\theta
\right],
\qquad
\theta = \frac{1}{4|\omega|\tau}.
Valid for all real :math:`\omega`. The function is even in :math:`\omega`,
non-negative, and satisfies:
- :math:`\tilde{\phi}(0) = 2\tau/\pi` (mean relaxation time)
- :math:`\tilde{\phi}(\omega) \propto \omega^{-1/2}` as :math:`\omega\tau \to \infty`
Note the elegant structure: :math:`C` and :math:`S` are both evaluated at
:math:`\theta`, the same argument as :math:`\sin\theta` and :math:`\cos\theta`.
**Fresnel Integral Convention**
The Fresnel integrals are defined here as (analog to ref 8 in [deGennes19]_)
.. math::
C(x) = \sqrt{\frac{2}{\pi}}\int_0^{\sqrt{x}} \cos(t^2)\,\mathrm{d}t,
\qquad
S(x) = \sqrt{\frac{2}{\pi}}\int_0^{\sqrt{x}} \sin(t^2)\,\mathrm{d}t,
with common limit
.. math::
C(\infty) = S(\infty) = \tfrac{1}{2}.
Relation to the SciPy convention :func:`scipy.special.fresnel`
(denoted :math:`C_s`, :math:`S_s`, which use :math:`\cos(\pi t^2/2)`):
.. math::
C(x) = C_s\!\left(\sqrt{\tfrac{2x}{\pi}}\right),
\qquad
S(x) = S_s\!\left(\sqrt{\tfrac{2x}{\pi}}\right).
The **Kohlrausch--Williams--Watts** (KWW) stretched exponential is
.. math::
\phi(t) = \exp\!\left[-\left(\frac{t}{\tau}\right)^{\!\beta}\right],
\qquad t \geq 0,
where :math:`\tau > 0` is the characteristic time and :math:`\beta = 1/2`.
The mean relaxation time (first relaxation time) is
.. math::
\langle\tau\rangle = \int_0^\infty \phi(t)\,\mathrm{d}t
= \frac{\tau}{\beta}\,\Gamma\!\left(\frac{1}{\beta}\right)
\;\xrightarrow{\;\beta=1/2\;}\; 2\tau.
The cosine Fourier transform convention used here is
.. math::
\tilde{\phi}(\omega)
= \int_0^\infty \cos(\omega t)\,\phi(t)\,\mathrm{d}t.
.. collapse:: Derivation
Step 1 — Substitution :math:`s=\sqrt{t}`
Set :math:`s = \sqrt{t}`, so :math:`t = s^2`, :math:`\mathrm{d}t = 2s\,\mathrm{d}s`
and :math:`\phi(t) = e^{-s/\sqrt{\tau}}`:
.. math::
\tilde{\phi}(\omega)
= 2\int_0^\infty s\,\cos(\omega s^2)\,e^{-s/\sqrt{\tau}}\,\mathrm{d}s.
Step 2 — Integration by Parts
Using :math:`\tfrac{\mathrm{d}}{\mathrm{d}s}\sin(\omega s^2) = 2\omega s\cos(\omega s^2)`,
the boundary term vanishes and
.. math::
\tilde{\phi}(\omega)
= \frac{1}{\omega\sqrt{\tau}}
\underbrace{
\int_0^\infty \sin(\omega s^2)\,e^{-s/\sqrt{\tau}}\,\mathrm{d}s
}_{K(\omega,\,\tau)}.
Step 3 — Reduction to a Laplace–Fresnel Integral
Substitute :math:`t = s\sqrt{\omega}`, so :math:`\omega s^2 = t^2` and
:math:`s/\sqrt{\tau} = \beta t` with
.. math::
\beta = \frac{1}{\sqrt{\omega\tau}}.
Then :math:`K = (1/\sqrt{\omega})\,\mathrm{Im}\,G(\beta)`, giving
.. math::
:label: result_G
\tilde{\phi}(\omega)
= \frac{\mathrm{Im}\,G(\beta)}{\omega^{3/2}\sqrt{\tau}},
\qquad
G(\beta) = \int_0^\infty e^{it^2 - \beta t}\,\mathrm{d}t.
Step 4 — Contour Evaluation of :math:`G(\beta)`
Complete the square in the exponent,
.. math::
it^2 - \beta t
= i\!\left(t + \tfrac{i\beta}{2}\right)^{\!2} + \tfrac{i\beta^2}{4},
and substitute :math:`v = t + i\beta/2`. The path of :math:`v` is a
horizontal line at height :math:`\mathrm{Im}(v) = \beta/2`. Since
:math:`|e^{iv^2}| = e^{-2\,\mathrm{Re}(v)\,\mathrm{Im}(v)} \to 0` as
:math:`\mathrm{Re}(v) \to \infty`, Cauchy's theorem allows deformation to
the real axis:
.. math::
G(\beta)
= e^{i\beta^2/4}
\left[
\int_0^\infty e^{iv^2}\,\mathrm{d}v
- \int_0^{i\beta/2} e^{iv^2}\,\mathrm{d}v
\right].
The first integral equals :math:`C(\infty) + iS(\infty) = (1+i)/2`.
For the imaginary-axis segment, parametrise :math:`v = iy` with
:math:`y: 0 \to \beta/2`, :math:`\mathrm{d}v = i\,\mathrm{d}y`:
.. math::
\int_0^{i\beta/2} e^{iv^2}\,\mathrm{d}v
= i\int_0^{\beta/2} e^{-iy^2}\,\mathrm{d}y
= S\!\left(\tfrac{\beta}{2}\right) + i\,C\!\left(\tfrac{\beta}{2}\right).
Since :math:`\beta/2 = \sqrt{\theta}` where :math:`\theta \equiv \beta^2/4`,
the Fresnel integrals with argument :math:`\beta/2` are related to those
with argument :math:`\theta` by
.. math::
\int_0^{\beta/2} \cos(t^2)\,\mathrm{d}t = \sqrt{\tfrac{\pi}{2}}\,C(\theta),
\qquad
\int_0^{\beta/2} \sin(t^2)\,\mathrm{d}t = \sqrt{\tfrac{\pi}{2}}\,S(\theta).
Assembling all pieces:
.. math::
G(\beta)
= e^{i\theta}
\left[
\left(\tfrac{1}{2} - S(\theta)\right)
+ i\left(\tfrac{1}{2} - C(\theta)\right)
\right].
**Final Result** in not normalized form
Taking :math:`\mathrm{Im}\,G(\beta)` and inserting into :eq:`result_G`:
.. math::
\tilde{\phi}(\omega)
= \sqrt{\frac{\pi}{2\omega^3\tau}}
\left[
\left(\tfrac{1}{2} - S(\theta)\right)\sin\theta
+
\left(\tfrac{1}{2} - C(\theta)\right)\cos\theta
\right],
\qquad
\theta = \frac{1}{4|\omega|\tau}.
Valid for all real :math:`\omega`. The function is even in :math:`\omega`,
non-negative, and satisfies:
- :math:`\tilde{\phi}(0) = 2\tau` (mean relaxation time)
- :math:`\tilde{\phi}(\omega) \propto \omega^{-1/2}` as :math:`\omega\tau \to \infty`
Note the elegant structure: :math:`C` and :math:`S` are both evaluated at
:math:`\theta`, the same argument as :math:`\sin\theta` and :math:`\cos\theta`.
**Asymptotic Limits**
**Low frequency**: :math:`\omega\tau \ll 1`
As :math:`\omega \to 0`, :math:`\theta \to \infty` and
:math:`C(\theta), S(\theta) \to 1/2`, so both bracket terms vanish.
The prefactor :math:`\omega^{-3/2}` is compensated and one recovers
.. math::
\tilde{\phi}(\omega) \;\xrightarrow{\;\omega\tau \to 0\;}\; 2\tau.
**High frequency**: :math:`\omega\tau \gg 1`
As :math:`\omega \to \infty`, :math:`\theta \to 0`,
:math:`C(\theta) \approx S(\theta) \approx 0`,
:math:`\cos\theta \to 1`, :math:`\sin\theta \to 0`, so
.. math::
\tilde{\phi}(\omega)
\;\xrightarrow{\;\omega\tau \to \infty\;}\;
\frac{1}{2}\sqrt{\frac{\pi\tau}{\omega}}
\;\propto\; \omega^{-1/2}.
This :math:`\omega^{-1/2}` power law is the hallmark of Rouse-mode dynamics
and arises naturally because :math:`\beta = 1/2` is produced by Fresnel-type
integrals over Gaussian chain modes (de Gennes, 1979).
References
-----------
.. [KWW1854] R. Kohlrausch,
*Theorie des elektrischen Rückstandes in der Leidener Flasche*,
Ann. Phys. Chem. **91**, 179 (1854).
.. [WW1970] G. Williams and D. C. Watts,
*Non-symmetrical dielectric relaxation behaviour arising from a simple
empirical decay function*,
Trans. Faraday Soc. **66**, 80 (1970).
.. [deGennes1979] P.-G. de Gennes,
*Scaling Concepts in Polymer Physics*,
Cornell University Press, Ithaca (1979).
.. [AS1972] M. Abramowitz and I. A. Stegun,
*Handbook of Mathematical Functions*,
Dover, New York (1972), Chap. 7.
.. [deGennes19] Quasi-Elastic Scattering of Neutrons by Dilute Polymer Solutions: I. Free-Draining Limit.
de Gennes, P.-G.
Physics Physique Fizika 1967, 3 (1), 37–45. https://doi.org/10.1103/PhysicsPhysiqueFizika.3.37.
"""
w = np.atleast_1d(np.asarray(w, dtype=float))
Y = np.zeros_like(w)
nonzero = np.abs(w) > 1e-10
aw = np.abs(w[nonzero])
tau = 1 / gamma
theta = 1.0 / (4.0 * aw * tau)
S, C = special.fresnel(np.sqrt(2 * theta / np.pi))
Y[nonzero] = (np.sqrt(np.pi / 2) / (aw**1.5 * np.sqrt(tau))
* ((0.5 - S) * np.sin(theta)
+ (0.5 - C) * np.cos(theta)))
Y[~nonzero] = 2.0 * tau
# with normalisation to integral = 1
result = dA(np.c_[w, Y / np.pi ].T)
result.setColumnIndex(iey=None)
result.tau = tau
result.columnname = 'w;kww'
result.modelname = inspect.currentframe().f_code.co_name
return result
[docs]
def stretchedExp_w(w, gamma, beta, tfactor=7):
r"""
KWW in frequency domain as Fourier tranform of stretched exponential.
.. math:: I(\omega) = FFT(I(t)) = FFT(e^{-(t\gamma)^\beta})
\text{ with normalization} \int_{-\infty}^{\infty} I(\omega) = 1
Parameters
----------
w : array
Frequencies in 1/ns
gamma : float
Relaxation rate in units 1/[unit t]
beta : float
Stretched exponent
tfactor : float
Power factor ``2**tfactor``to expand time range and time steps in FFT to improve accuracy.
Returns
-------
dataArray
Examples
--------
::
import jscatter as js
import numpy as np
w = np.r_[-80:80:0.1]
p = js.grace()
g=2
sqw1 = js.dynamic.stretchedExp_w(w, gamma=g, beta=1)
sqwL = js.dynamic.lorentz_w(w, hwhm = g)
maxerr = (np.abs(sqw1.Y-sqwL.Y)/sqwL.Y).max()
p.plot(sqw1,sy=[1,0.4,1],le=fr'stretched beta=1 max rel err :{maxerr:.4g}')
p.plot(sqwL,sy=0,li=[1,1,11],le='Lorentz analytic')
sqw05 = js.dynamic.stretchedExp_w(w, gamma=g, beta=0.5,tfactor=23)
sqw05a = js.dynamic.stretchedExp_beta05_w(w,g)
maxerr = (np.abs(sqw05.Y-sqw05a.Y)/sqw05a.Y).max()
p.plot(sqw05,sy=[1,0.4,2],le=fr'stretched beta=0.5 max rel err :{maxerr:.4g}')
p.plot(sqw05a,sy=0,li=[1,1,5],le='stretched beta=0.5 analytic')
p.yaxis(label='S(w)',scale='log',min=1e-4,max=1)
p.xaxis(label='w / 1/ns',scale='norm',min=-80,max=80)
p.legend(x=15,y=sqw05.Y.max(),charsize=0.7)
p.title('Comparison FFT(stretched exp.) to analytic result')
# p.save(js.examples.imagepath+'/stretchedExp_w.png',size=(1.5,1.5),dpi=300)
.. image:: ../../examples/images/stretchedExp_w.png
:align: center
:width: 60 %
:alt: ORZ eigenvalue and more
"""
tfactor = int(tfactor)
w = np.atleast_1d(w)
# determine the times and differences dt
dt = 1 / np.max(np.abs(w)) / tfactor
nn = 2**(int(np.log2(w.shape[0] * tfactor)))
tt = np.r_[0:nn] * dt
Y = np.exp(-(tt * gamma) ** beta)
# make it symmetric zero only once
RY = np.r_[Y[:0:-1], Y]
# do rfft from -N to N
# using spectrum from -N,N the shift theorem says we get a
# exp[-j*2*pi*f*N/2] phase leading to alternating sign => use the absolute value
# the 2*pi prefactor comes from the scipy/numpy FFT convention
# which uses exp(-2pi*j...)[technical convention] and not exp(-j)[physics convention]
wn = 2 * pi * np.fft.rfftfreq(2 * nn - 1, dt) # frequencies
wY = dt * np.abs(np.fft.rfft(RY).real) / (2 * pi) # fft
# now try to average or interpolate for needed w values
wn = np.r_[-wn[:0:-1], wn]
wY = np.r_[wY[:0:-1], wY]
# integral = scipy.integrate.simpson(y=wY, x=wn)
#---------------------
result = dA(np.c_[w, np.interp(w,wn,wY)].T)
# result.Sq = integral
result.setColumnIndex(iey=None)
result.columnname = 'w;Iqw;I'
result.gamma = gamma
result.beta = beta
result.modelname = inspect.currentframe().f_code.co_name
return result
_erfi = special.erfi
_G = special.gamma
_h1f1 = special.hyp1f1
_erf = special.erf
_Gi = special.gammainc
def _ln(w, n):
return n / pi / (n * n + w * w)
[docs]
def lorentz_w(w, hwhm):
r"""
Normalized Lorentz function.
.. math :: Ln(w) = \frac{\gamma}{\pi(w^2+\gamma^2)}
Same as KWW in frequency domain as Fourier tranform of stretched exponential with beta=1.
Parameters
----------
w : array
Frequency
hwhm : float
Half width half maximum :math:`\gamma` or relaxation rate.
Returns
-------
dataArray
Notes
-----
Examples
--------
::
import jscatter as js
import numpy as np
w = np.r_[-100:100:0.1]
p = js.grace(2,2)
for hw in np.r_[1:32:5]:
iw = js.dynamic.lorentz_w(w,hw)
p.plot(iw,le=f'hw={hw}')
p.yaxis(scale='l',label=r'S(\xw\f{}) / a.u.',min=1e-5,max=1)
p.xaxis(label=r'\xw\f{} / ns\S-1',min=-100,max=100)
p.legend(x=30,y=1,charsize=0.8)
# p.save(js.examples.imagepath+'/lorentz_w.jpg', size=(2,2))
.. image:: ../../examples/images/lorentz_w.jpg
:align: center
:width: 50 %
:alt: dynamic_t2f_examples.
"""
result = dA(np.c_[w, _ln(w, hwhm)].T)
result.modelname = inspect.currentframe().f_code.co_name
result.setColumnIndex(iey=None)
result.columnname = 'w;Ln'
result.hwhm = hwhm
return result
[docs]
def diffusionHarmonicPotential_w(w, q, tau, rmsd, ndim=3, nmax='auto'):
r"""
Diffusion in a harmonic potential for dimension 1,2,3 (isotropic averaged), dynamic structure factor in w domain.
An approach worked out by Volino et al.[1]_ assuming Gaussian confinement and leads to a more efficient
formulation by replacing the expression for diffusion in a sphere with a simpler expression pertaining
to a soft confinement in harmonic potential.
:math:`D_t = \langle u_x^2 \rangle / \tau_0` see equ. 32 in [1]_ .
Parameters
----------
w : array
Frequencies in 1/ns
q : float
Wavevector in nm**-1
tau : float
Mean correlation time :math:`\tau_0`. In units ns.
rmsd : float
Root mean square displacement :math:`\langle u_x^2 \rangle^{1/2}` (width) of the Gaussian in units nm.
ndim : 1,2,3, default=3
Dimensionality of the potential.
nmax : int,'auto'
Order of expansion.
'auto' -> nmax = min(max(int(6*q * q * u2),30),1000)
Returns
-------
dataArray
Notes
-----
Volino et al.[1]_ compared the behavior of this approach to the well known expression for diffusion in a sphere.
Even if the details differ, the salient features of both models match if the radius R**2 ≃ 5*u0**2 and
the diffusion constant inside the sphere relates to the relaxation time of particle correlation t0= ⟨u**2⟩/Ds
towards the Gaussian with width u0=⟨u**2⟩**0.5.
.. math:: I_s(Q_x,\omega) = A_0(Q) + \sum_n^{\infty} A_n(Q) L_n(\omega)
\; with \; L_n(\omega) = \frac{\tau_0 n}{\pi (n^2+ \omega^2\tau_0^2)}
ndim=3
Here we use the Fourier transform of equ 23 with equ. 27a+b in [1]_.
For order n>30 the Stirling approximation for n! in equ 27b of [1]_ is used.
.. math:: A_0(Q) = e^{-Q^2\langle u^2_x \rangle}
.. math:: A_n(Q,\omega) = e^{-Q^2\langle u^2_x \rangle} \frac{(Q^2\langle u^2_x \rangle)^n}{n!}
ndim=2
Here we use the Fourier transform of equ 23 with equ. 28a+b in [1]_.
.. math:: A_0(Q) = \frac{\sqrt{\pi} e^{-Q^2\langle u^2_x \rangle}}{2}
\frac{erfi(\sqrt{Q^2\langle u^2_x \rangle})}{\sqrt{Q^2\langle u^2_x \rangle}}
.. math:: A_n(Q,\omega) = \frac{\sqrt{\pi} (Q^2\langle u^2_x \rangle)^n}{2}
\frac{F_{1,1}(1+n;3/2+n;-Q^2\langle u^2_x \rangle)}{\Gamma(3/2+n)}
with :math:`F_{1,1}(a,b,z)` Kummer confluent hypergeometric function, Gamma function :math:`\Gamma`
and *erfi* is the imaginary error function *erf(iz)/i*
ndim=1
The equation given by Volino (29a+b in [1]_) seems to be wrong as a comparison with the Fourier transform and
the other dimensions shows.
Use the model from time domain and use FFT as shown in the example.
For experts: To test this remove a flag in the source code and compare.
Examples
--------
::
import jscatter as js
import numpy as np
t2f = js.dynamic.time2frequencyFF
dHP = js.dynamic.diffusionHarmonicPotential
w = np.r_[-100:100]
ql = np.r_[1,3,6,9,12,15]
iqt3 = js.dL([js.dynamic.diffusionHarmonicPotential_w(w=w,q=q,tau=0.14,rmsd=0.34,ndim=3) for q in ql])
iqt2 = js.dL([js.dynamic.diffusionHarmonicPotential_w(w=w,q=q,tau=0.14,rmsd=0.34,ndim=2) for q in ql])
# as ndim=1 is a wrong solution use this instead
# To move spectral leakage out of our window we increase w and interpolate.
# The needed factor (here 23) depends on the quality of your data and background contribution.
# You may test it using ndim=2 in this example.
iqt1 = js.dL([t2f(dHP,'elastic',w=w*23,q=q, rmsd=0.34, tau=0.14 ,ndim=1).interpolate(w) for q in ql])
p=js.grace(1,1)
p.multi(2,3)
p[1].title('diffusionHarmonicPotential for ndim= 1,2,3')
for i,(i3,i2,i1) in enumerate(zip(iqt3,iqt2,iqt1)):
p[i].plot(i3,li=1,sy=0,le='$wavevector nm\S-1')
p[i].plot(i2,li=2,sy=0)
p[i].plot(i1,li=4,sy=0)
p[i].yaxis(scale='log')
if i in [1,2,4,5]:p[i].yaxis(ticklabel=0)
p[i].legend(x=5,y=1, charsize=0.7)
p[0].xaxis(label='')
p[0].yaxis(label=r'S(\xw)')
p[3].yaxis(label=r'S(\xw)')
# p.save(js.examples.imagepath+'/diffusionHarmonicPotential.jpg', size=(2,2))
.. image:: ../../examples/images/diffusionHarmonicPotential.jpg
:align: center
:width: 50 %
:alt: dynamic_t2f_examples.
References
----------
.. [1] Gaussian model for localized translational motion: Application to incoherent neutron scattering.
Volino, F., Perrin, J. C. & Lyonnard, S. J. Phys. Chem. B 110, 11217–11223 (2006).
"""
w = np.array(w, float)
u2 = rmsd ** 2
if not isinstance(nmax, numbers.Integral):
nmax = min(max(int(6 * q * q * u2), 30), 1000)
Ln = lambda w, t0, n: t0 / pi * n / (n * n + w * w * t0 * t0) # equ 25a
if ndim == 3:
# 3D case
A0 = lambda q: np.exp(-q * q * u2) # EISF equ 27a
def An(q, n):
s = (n < 30) # select not to large n and use for the other the Stirling equation
An = np.r_[
(q * q * u2) ** n[s] / special.factorial(n[s]), (q * q * u2 / n[~s] * np.e) ** n[~s] / (
2 * pi * n[~s]) ** 0.5]
An *= np.exp(-q * q * u2)
return An
n = np.r_[:nmax] + 1
an = An(q, n)
sel = np.isfinite(an) # remove An with inf or nan
Iqw = (an[sel, None] * Ln(w, tau, n[sel, None])).sum(axis=0) # equ 23 after ft
Iqw[np.abs(w) < 1e-8] += A0(q)
elif ndim == 2:
# 2D case
A0 = lambda q: pi ** 0.5 / 2. * np.exp(-q * q * u2) * _erfi((q * q * u2) ** 0.5) / (
q * q * u2) ** 0.5 # EISF equ 28a
An = lambda q, n: pi ** 0.5 / 2. * (q * q * u2) ** n * _h1f1(1 + n, 1.5 + n, -q * q * u2) / _G(
1.5 + n) # equ 28b
n = np.r_[:nmax] + 1
Iqw = (An(q, n)[:, None] * Ln(w, tau, n[:, None])).sum(axis=0) # equ 23 after ft
Iqw[np.abs(w) < 1e-8] += A0(q)
elif ndim == 1 and False:
print(' THis seems to be wrong as given in the paper')
# 1D case
A0 = lambda q: pi ** 0.5 / 2. * _erf((q * q * u2) ** 0.5) / (q * q * u2) ** 0.5 # EISF equ 29a
An = lambda q, n: (_G(0.5 + n) - _Gi(0.5 + n, q * q * u2)) / (2 * (q * q * u2) ** 0.5 * _G(1 + n)) # equ 29b
n = np.r_[:nmax] + 1
an = An(q, n)
sel = np.isfinite(an) # remove An with inf or nan
Iqw = (an[sel, None] * Ln(w, tau, n[sel, None])).sum(axis=0) # equ 23 after ft
Iqw[np.abs(w) < 1e-8] += A0(q)
else:
raise Exception('ndim should be one of 2 or 3; for 1 use fourier tranform from time domain, see doc.')
result = dA(np.c_[w, Iqw].T)
result.modelname = inspect.currentframe().f_code.co_name
result.setColumnIndex(iey=None)
result.columnname = 'w;Iqw'
result.u0 = rmsd
result.dimension = ndim
result.wavevector = q
result.meancorrelationtime = tau
result.gaussWidth = rmsd
result.nmax = nmax
result.Ds = rmsd ** 2 / tau
return result
#: First 99 coefficients from Volino for diffusionInSphere_w
# VolinoCoefficient=np.loadtxt(os.path.join(_path_,'data','VolinoCoefficients.dat')) # numpy cannot load because of utf8
with open(os.path.join(_path_, '../data', 'VolinoCoefficients.dat')) as f: VolinoC = f.readlines()
VolinoCoefficient = np.array([line.strip().split() for line in VolinoC if line[0] != '#'], dtype=float)
[docs]
def diffusionInSphere_w(w, q, D, R):
r"""
Diffusion inside of a sphere; dynamic structure factor in w domain.
Parameters
----------
w : array
Frequencies in 1/ns
q : float
Wavevector in nm**-1
D : float
Diffusion coefficient in units nm**2/ns
R : float
Radius of the sphere in units nm.
Returns
-------
dataArray
Notes
-----
Here we use equ. 33 in [1]_
.. math:: S(q,\omega) = A_0^0(q) \delta(\omega) + \frac{1}{\pi}
\sum_{l,n\ne 0,0}(2l+1)A_n^l(q) \frac{(x_n^l)^2D/a^2}{[(x_n^l)^2D/a^2]^2 + \omega^2}
with :math:`x_n^l` as the first 99 solutions of equ 27 a+b as given in [1]_ and
.. math:: A_0^0(q) = \big[ \frac{3j_1(qa)}{qa} \big]^2 , \; (l,n) = (0,0)
.. math:: A_n^l(q) &= \frac{6(x_n^l)^2}{(x_n^l)^2-l(l+1)}
\big[\frac{qaj_{l+1}(qa)-lj_l(qa)}{(qa)^2-(x_n^l)^2}\big]^2 \; for \; qa\ne x_n^l
&= \frac{3}{2}j_l^2(x_n^l) \frac{(x_n^l)^2-l(l+1)}{(x_n^l)^2} \; for \; qa = x_n^l
This is valid for qR<20 with accuracy of ~0.001 as given in [1]_.
If we look at a comparison with free diffusion the valid range seems to be smaller.
A comparison of diffusion in different restricted geometry is show in example
:ref:`A comparison of different dynamic models in frequency domain`.
Examples
--------
::
import jscatter as js
import numpy as np
w=np.r_[-100:100]
ql=np.r_[1:14.1:1.3]
p=js.grace(1,1)
iqw=js.dL([js.dynamic.diffusionInSphere_w(w=w,q=q,D=0.14,R=0.2) for q in ql])
p.plot(iqw)
p.yaxis(label=r'S(\xw)',scale='l')
p.xaxis(label=r'\xw\f{} / ns\S-1')
# p.save(js.examples.imagepath+'/diffusionInSphere_w.jpg', size=(2,2))
.. image:: ../../examples/images/diffusionInSphere_w.jpg
:align: center
:width: 50 %
:alt: dynamic_t2f_examples.
References
----------
.. [1] Neutron incoherent scattering law for diffusion in a potential of spherical symmetry:
general formalism and application to diffusion inside a sphere.
Volino, F. & Dianoux, A. J., Mol. Phys. 41, 271–279 (1980).
https://doi.org/10.1080/00268978000102761
"""
nmax = 99
qR = q * R
x = VolinoCoefficient[1:nmax, 0] # x_n_l
x2 = x ** 2
l = VolinoCoefficient[1:nmax, 1].astype(int)
# n = VolinoCoefficient[1:50, 2].astype(int)
w = np.array(w, float)
Ln = lambda w, g: g / (g * g + w * w)
A0 = lambda qa: (3 * spjn(1, qa) / qa) ** 2
def Anl(qa):
# equ 31 a+b in [1]_
res = np.zeros_like(x)
s = (x == qa)
if np.any(s):
res[s] = 1.5 * spjn(l[s], x[s]) ** 2 * (x2[s] - l[s] * (l[s] + 1)) / x2[s]
if np.any(~s):
s = ~s # not s
res[s] = 6 * x2[s] / (x2[s] - l[s] * (l[s] + 1)) * (
(qa * spjn(l[s] + 1, qa) - l[s] * spjn(l[s], qa)) / (qa ** 2 - x2[s])) ** 2
return res
Iqw = 1 / pi * (((2 * l + 1) * Anl(qR))[:, None] * Ln(w, x2[:, None] * D / R ** 2)).sum(axis=0) # equ 33
Iqw[np.abs(w) < 1e-8] += A0(q)
result = dA(np.c_[w, Iqw].T)
result.modelname = inspect.currentframe().f_code.co_name
result.setColumnIndex(iey=None)
result.columnname = 'w;Iqw'
result.radius = R
result.wavevector = q
result.diffusion = D
return result
[docs]
def rotDiffusion_w(w, q, cloud, Dr, lmax='auto'):
r"""
Rotational diffusion of an object (dummy atoms); dynamic structure factor in w domain.
A cloud of dummy atoms can be used for coarse graining of a non-spherical object e.g. for amino acids in proteins.
On the other hand its just a way to integrate over an object e.g. a sphere or ellipsoid.
We use [2]_ for an objekt of arbitrary shape modified for incoherent scattering.
Parameters
----------
w : array
Frequencies in 1/ns
q : float
Wavevector in units 1/nm
cloud : array Nx3, Nx4 or Nx5 or float
- A cloud of N dummy atoms with positions cloud[:3] that describe an object.
- If given, cloud[3] is the incoherent scattering length :math:`b_{inc}` otherwise its equal 1.
- If given, cloud[4] is the coherent scattering length otherwise its equal 1.
- If cloud is single float the value is used as radius of a sphere with 10x10x10 grid.
Dr : float
Rotational diffusion constant in units 1/ns.
lmax : int
Maximum order of spherical bessel function.
'auto' -> lmax > 2pi*r.max()*q/6.
Returns
-------
dataArray
Columns [w; Iqwinc; Iqwcoh]
Input parameters as attributes.
Notes
-----
See :py:func:`~.timedomain.transRotDiffusion` for more details.
The Fourier transform of the *exp* function is a Lorentzian so the *exp* should be changed to a Lorentzian.
Examples
--------
::
import jscatter as js
import numpy as np
R=2;NN=10
Drot=js.formel.Drot(R)
ql=np.r_[0.5,2:18.:3]
w=np.r_[-100:100:0.1]
grid=js.ff.superball(ql,R,p=1,nGrid=NN,returngrid=True)
p=js.grace()
iqwR1=js.dL([js.dynamic.rotDiffusion_w(w,q,grid.XYZ,Drot) for q in ql])
p.plot(iqwR1,le=f'NN={NN:.0f} q=$wavevector nm\S-1')
p.yaxis(scale='l',label=r'S(\xw\f{}) / a.u.',min=1e-4,max=1e4)
p.xaxis(label=r'\xw\f{} / ns\S-1',min=-100,max=100)
p.legend(x=30,y=9000,charsize=0.8)
# p.save(js.examples.imagepath+'/transRotDiffusion_w.jpg', size=(2,2))
.. image:: ../../examples/images/transRotDiffusion_w.jpg
:align: center
:width: 50 %
:alt: dynamic_t2f_examples.
References
----------
.. [1] Incoherent scattering law for neutron quasi-elastic scattering in liquid crystals.
Dianoux, A., Volino, F. & Hervet, H. Mol. Phys. 30, 37–41 (1975).
.. [2] Effect of rotational diffusion on quasielastic light scattering from fractal colloid aggregates.
Lindsay, H., Klein, R., Weitz, D., Lin, M. & Meakin, P. Phys. Rev. A 38, 2614–2626 (1988).
"""
#: Lorentzian
Ln = lambda w, g: g / (g * g + w * w) / pi
if isinstance(cloud, numbers.Number):
R = cloud
NN = 10
grid = np.mgrid[-R:R:1j * NN, -R:R:1j * NN, -R:R:1j * NN].reshape(3, -1).T
inside = lambda xyz, R: la.norm(grid, axis=1) < R
cloud = grid[inside(grid, R)]
if cloud.shape[1] == 5:
# last columns are incoherent and coherent scattering length
blinc = cloud[:, 3]
blcoh = cloud[:, 4]
cloud = cloud[:, :3]
elif cloud.shape[1] == 4:
# last column is scattering length
blinc = cloud[:, 3]
blcoh = np.ones(cloud.shape[0])
cloud = cloud[:, :3]
else:
blinc = np.ones(cloud.shape[0])
blcoh = blinc
w = np.array(w, float)
bi2 = blinc ** 2
r, p, t = formel.xyz2rphitheta(cloud).T
pp = p[:, None]
tt = t[:, None]
qr = q * r
if not isinstance(lmax, numbers.Integral):
# lmax = pi * r.max() * q / 6. # a la CRYSON (SANS/SAXS)
# we need a factor of 2 more compared to CRYSON for Q>10 nm**-1
lmax = min(max(2 * int(pi * qr.max() / 6. * 2), 7), 100)
# We calc here the field autocorrelation function as in equ 24
# Fourier transform of the exp result in lorentz function
# incoherent with i=j -> Sum_m(Ylm) leads to (2l+1)/4pi
bjlylminc = [(bi2 * spjn(l, qr) ** 2 * (2 * l + 1)).sum() for l in np.r_[:lmax + 1]]
# add Lorentzian
Iqwinc = np.c_[[bjlylminc[l].real * Ln(w, l * (l + 1) * Dr) for l in np.r_[:lmax + 1]]].sum(axis=0)
Iq_inc = np.sum(bjlylminc).real
# coh is sum over i then squared and sum over m see Lindsay equ 19
bjlylmcoh = [4 * pi * np.sum(np.abs((blcoh * spjn(l, qr) * Ylm(l, np.r_[-l:l + 1], pp, tt).T).sum(axis=1)) ** 2)
for l in np.r_[:lmax + 1]]
Iqwcoh = np.c_[[bjlylmcoh[l].real * Ln(w, l * (l + 1) * Dr) for l in np.r_[:lmax + 1]]].sum(axis=0)
Iq_coh = np.sum(bjlylmcoh).real
result = dA(np.c_[w, Iqwinc, Iqwcoh].T)
result.modelname = inspect.currentframe().f_code.co_name
result.setColumnIndex(iey=None)
result.columnname = 'w; Iqwinc; Iqwcoh'
result.radiusOfGyration = np.sum(r ** 2) ** 0.5
result.Iq_coh = Iq_coh
result.Iq_inc = Iq_inc
result.wavevector = q
result.rotDiffusion = Dr
result.lmax = lmax
return result
[docs]
def nSiteJumpDiffusion_w(w, q, N, t0, r0):
r"""
Random walk among N equidistant sites (isotropic averaged); dynamic structure factor in w domain.
E.g. for CH3 group rotational jump diffusion over 3 sites.
Parameters
----------
w : array
Frequencies in 1/ns
q: float
Wavevector in units 1/nm
N : int
Number of jump sites, jump angle 2pi/N
r0 : float
Distance of sites from center of rotation.
For CH3 e.g.0.12 nm.
t0 : float
Rotational correlation time.
Returns
-------
dataArray
Notes
-----
Equ. 24 [1]_ :
.. math:: S_{inc}^{rot}(Q,\omega) = B_0(Qa)\delta(\omega) + \frac{1}{\pi} \sum_{n=1}^{N-1} B_n(Qa)
\frac{\tau_n}{1+(\omega\tau_n)^2}
with :math:`\tau_1=\frac{\tau}{1-cos(2\pi/N)}` , :math:`\tau_n=\tau_1\frac{sin^2(\pi/N)}{sin^2(n\pi/N)}`
.. math:: B_n(Qa) = \frac{1}{N} \sum_{p=1}^{N} j_0 \Big( 2Qa sin(\frac{\pi p}{N}) \Big) cos(n\frac{2\pi p}{N})
Examples
--------
::
import jscatter as js
import numpy as np
w=np.r_[-100:100:0.1]
ql=np.r_[1:14.1:1.3]
p=js.grace()
iqw=js.dL([js.dynamic.nSiteJumpDiffusion_w(w=w,q=q,N=3,t0=0.01,r0=0.12) for q in ql])
p.plot(iqw)
p.yaxis(scale='l',label=r'S(\xw\f{}) / a.u.',min=1e-6,max=1)
p.xaxis(label=r'\xw\f{} / ns\S-1',min=-100,max=100)
# p.save(js.examples.imagepath+'/nSiteJumpDiffusion_w.jpg', size=(2,2))
.. image:: ../../examples/images/nSiteJumpDiffusion_w.jpg
:align: center
:width: 50 %
:alt: dynamic_t2f_examples.
References
----------
.. [1] Incoherent scattering law for neutron quasi-elastic scattering in liquid crystals.
Dianoux, A., Volino, F. & Hervet, H., Mol. Phys. 30, 37–41 (1975).
https://doi.org/10.1080/00268977500102721
"""
w = np.array(w, float)
#: Lorentzian
Ln = lambda w, tn: tn / (1 + (w * tn) ** 2) / pi
def Bn(qa, n):
return np.sum([spjn(0, 2 * qa * np.sin(pi * p / N)) * np.cos(n * 2 * pi * p / N) for p in np.r_[:N] + 1]) / N
B0 = np.sum([spjn(0, 2 * q * r0 * np.sin(pi * p / N)) for p in np.r_[:N] + 1]) / N
t1 = t0 / (1 - np.cos(2 * pi / N))
tn = lambda n: t1 * np.sin(pi / N) ** 2 / np.sin(n * pi / N) ** 2
Iqw = np.c_[[Bn(q * r0, n) * Ln(w, tn(n)) for n in np.r_[1:N]]].sum(axis=0)
Iqw[np.abs(w) < 1e-8] += B0
result = dA(np.c_[w, Iqw].T)
result.modelname = inspect.currentframe().f_code.co_name
result.setColumnIndex(iey=None)
result.columnname = 'w;Iqw'
result.r0 = r0
result.wavevector = q
result.t0 = t0
result.N = N
return result
def _gauss(x, mean, sigma):
return np.exp(-0.5 * ((x - mean) / sigma) ** 2)
# noinspection PyIncorrectDocstring
[docs]
def resolution_w(w, s0=1, m0=0, s1=None, m1=None, s2=None, m2=None, s3=None, m3=None,
s4=None, m4=None, s5=None, m5=None,s6=None, m6=None, s7=None, m7=None,
a0=1, a1=1, a2=1, a3=1, a4=1, a5=1, a6=1, a7=1, bgr=0, resolution=None):
r"""
Resolution as multiple Gaussians for inelastic measurement as backscattering or time of
flight instrument in w domain.
Multiple Gaussians define the function to describe a resolution measurement.
Use only a common mi to account for a shift.
See :py:func:`~.dynamic.timedomain.resolution` for transform to time domain.
Parameters
----------
w : array
Frequencies
s0,s1,... : float
Sigmas of several Gaussian functions representing a resolution measurement.
The number of si not none determines the number of Gaussians.
m0, m1,.... : float, None
Means of the Gaussian functions representing a resolution measurement.
a0, a1,.... : float, None
Amplitudes of the Gaussian functions representing a resolution measurement.
bgr : float, default=0
Background
resolution : dataArray
Resolution with attributes sigmas, amps which are used instead of si, ai.
- If from t domain this represents the Fourier transform from w to t domain.
The means are NOT used from as these result only in a phase shift, instead m0..m5 are used.
- If from w domain the resolution is recalculated.
Returns
-------
dataArray :
- .means
- .amps
- .sigmas
Notes
-----
In a typical inelastic experiment the resolution is measured by e.g. a vanadium measurement (elastic scatterer).
This is described in `w` domain by a multi Gaussian function as in resw=resolution_w(w,...) with
amplitudes :math:`a_{iw}`, width sigma :math:`s_{iw}` and common mean :math:`m_w`.
To allow asymmetric resolutions as observed on some instruments we use mean :math:`m_{iw}`
resolution(t,resolution_w=resw) defines the Fourier transform of resolution_w using the same coefficients.
:math:`m_{it}` are set by default to 0 (if not explicit set) as :math:`m_{iw}` lead only to a phase shift.
It is easiest to shift w values in w domain as it corresponds to a shift of the elastic line.
The used Gaussians are normalized that they are a pair of Fourier transforms:
.. math:: R_t(t,m_i,s_i,a_i)=\sum_i a_i s_i e^{-\frac{1}{2}s_i^2 t^2} \Leftrightarrow
R_w(w,m_i,s_i,a_i)=\sum_i a_i e^{-\frac{1}{2}(\frac{w-m_i}{s_i})^2}
under the Fourier transform defined as
.. math:: F(f(t)) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} f(t) e^{-i\omega t} dt
.. math:: F(f(w)) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} f(\omega) e^{i\omega t} d\omega
Examples
--------
Transform from and to time domain
::
import jscatter as js
import numpy as np
# resw is a resolution in w domain maybe as a result from a fit to vanadium data
# resw contains all parameters
w=np.r_[-100:100:0.5]
resw=js.dynamic.resolution_w(w, s0=12, m0=0, a0=2)
w2=np.r_[0:50:0.2]
rest2=js.dynamic.resolution_w(w2,resolution=resw)
# representing the Fourier transform of to time domain
t=np.r_[0:1:0.01]
rest=js.dynamic.resolution(t,resolution=resw)
**Sequential fit in w domain to a measurement with real data.**
The data file is from the SPHERE instrument at MLZ Garching (usually not gzipped).
We overfitt the data as we just want a description of Gaussians with :math:`\chi^2 \approx 1`
::
import jscatter as js
import numpy as np
def readinx(filename,block, l2p):
# reusable function for inx reading
# neutron scatterers use strange data formats like .inx dependent on local contact
# in the beginning we find sometimes an index or the respective Q value or number of channels.
# here the blocks start with number of channels. Spaces (' 562 ') are important to find the starting line.
# this might dependent on instrument and local contact (they can do what you ask for)
data = js.dL(filename,
block=block, # this finds the starting of the blocks of all angles or Q
usecols=[1,2,3], # ignore the numbering at each line
lines2parameter=l2p) # catch the parameters at beginning
for da in data:
da.channels = da.line_1[0]
da.Q = da.line_3[0] # in 1/A
da.incident_energy = da.line_3[1] # in meV
da.temperature = da.line_3[3] # in K if given
return data
vanae = readinx(js.examples.datapath +'/Vana.inx.gz',block=' 562 ',l2p=[-1,-2,-3,-4])
# convert to 1/ns, we can select parts of the data
vanat = js.dynamic.convert_e2w(vanae[:],T=293,unit='μeV')
vana = js.dynamic.shiftAndBinning(vanat) # if you need the center w0
# low Q might be fitted with less Gaussians
# high Q with 6 Gaussians increasing width , decreasing amplitude
dm = 2
start = {'s0':0.3,'m0':0,'a0':max,'bgr':0.73,
's1':2,'m1':0,'a1':0.3*max,
's2':4,'m2':0,'a2':0.1*max,
's3':8,'m3':0,'a3':0.01*max,
's4':16,'m4':0,'a4':0.001*max,
's5':32,'m5':0,'a5':0.0001*max}
for van in vana:
max= van.Y.max()
van.setlimit(m0=[-dm,dm],m1=[-dm,dm],m2=[-dm,dm],m3=[-dm,dm],m4=[-dm,dm],m5=[-dm,dm])
van.setlimit(a0=[0],a1=[0],a2=[0],a3=[0],a4=[0],a5=[0])
van.setlimit(s0=[0],s1=[0],s2=[0],s3=[0],s4=[0],s5=[0])
van.makeErrPlot(yscale='log', fitlinecolor=11,title=f'Q={van.Q:.3f}')
van.fit(js.dynamic.resolution_w,start,{},{'w':'X'},
max_nfev=20000,method='lm',diff_step=0.001,workers=0)
van.lastfit.Q = van.Q # needed later
start = van.getFreepar # use last result as new start parameter
# save to recover later (that we dont need to repeat this)
vanalastfit = js.dL([v.lastfit for v in vana])
vanalastfit.save('Vana_fitted.inx.gz')
Recover the result to use for convolution in a later step
::
# recover by
vanalastfit = js.dL('Vana_fitted.inx.gz')
vanalastfit.getfromcomment('modelname')
# select a specific Q and recalculate; w can also be adopted to be different
vanQ = vanalastfit.filter(Q=1.417)[0]
res = js.dynamic.resolution_w(w=vanQ.X, resolution=vanQ)
# vana[6].showlastErrPlot(yscale='log',fitlinecolor=6,title=f'Q={vana[7].Q:.3f}')
# vana[6].savelastErrPlot(js.examples.imagepath+'/resolutionfit.jpg')
.. image:: ../../examples/images/resolutionfit.jpg
:align: center
:width: 50 %
:alt: worm
"""
if resolution is None:
means = [m0, m1, m2, m3, m4, m5, m6, m7]
sigmas = [s0, s1, s2, s3, s4, s5, s6, s7]
amps = [a0, a1, a2, a3, a4, a5, a6, a7]
else:
if resolution.modelname[-1] == 'w':
# resolution from w domain
means = resolution.means
sigmas = resolution.sigmas
amps = resolution.amps
else:
means = [0 if m is None else m for m in [m0, m1, m2, m3, m4, m5, m6, m7]]
sigmas = [1. / s if s is not None else s for s in resolution.sigmas]
amps = resolution.amps
w = np.atleast_1d(w)
if isinstance(resolution, str): # elastic
Y = elastic_w(w).Y
integral = 1
else:
# filter none numbers
sma = np.array([[s, m, a] for s, m, a in zip(sigmas, means, amps)
if (np.issubdtype(np.array([s, m, a]).dtype, np.number)) and (None not in [s, m, a])]).T
Y = np.sum(sma[2][:, None] * _gauss(x=w, mean=sma[1][:, None], sigma=sma[0][:, None]), axis=0)
integral = integrate.trapezoid(Y, w)
result = dA(np.c_[w, Y + bgr].T)
result.setColumnIndex(iey=None)
result.modelname = inspect.currentframe().f_code.co_name
result.columnname = 'w;Rw'
result.means = means
result.sigmas = sigmas
result.amps = amps
result.integral = integral
return result
[docs]
def doubleStretchedExp_w(w, Q, g1, beta1, amp1, resolution, g2=1, beta2=1, amp2=0, elastic=0, bgr=0, wmax=None):
r"""
Two stretched exponential with elastic contribution and resolution smearing.
A convenience function that allows to directly fit experimental data.
Remember that fixing beta can mimic different model like beta=1 is Lorentz function,
beta = 0.5 fits to subdiffusive Rouse like motions.
Each stretchedExp is normalized :math:`\int_{-\infty}^{\infty} \phi(\omega) = 1`
Parameters
----------
w : array
Frequencies 1/ns
Q : float
Wavevector in 1/nm for selecting the correct resolution.
g1,g2 : float
Relaxation rates in units 1/[unit t]
beta1,beta2 : float
Stretched exponent corresponding to g1,g2
amp1,amp2 : float
Amplitudes of the FF transformed double Exp.
elastic : float
Amplitude elastic contribution.
If the resolution is normalized it is directly the elastic contribution.
resolution : dataList
Measured or calculated resolution function with same Q as in the data to fit.
- Normalize the resolution before usage that the elastic amplitude a0 gets meaningful.
- In general it is assumed that the resolution and data w are equidistant.
- A fitted resolution can be used to smooth noise in a measured resolution.
See :py:func:`~.dynamic.frequencydomain.resolution_w` how to fit and use the '.lastfit' .
bgr : float
Background.
wmax : float, default None
Cutoff frequency for resolution. Not needed for fitted resolution.
Return
------
model : dataArray for Q
Notes
-----
We sum here two stretched exponentials with amplitudes and add an elastic contribution.
For stretched exp see :py:func:`stretchedExp_w` that are convoluted with the resolution cut at wmax.
As elastic contribution we use the measured resolution cut at wmax and multiplied by amplitude `elastic`
A constant bgr is added.
Examples
--------
::
import jscatter as js
import numpy as np
def readEMUdat(filename, wavelength=6.27084, temperature=0):
# Read data form EMU@ANSTO as direct output from Mantid (No .inx converison needed)
data = js.dL(filename, delimiter=',')
for da in data:
da.Q = np.round(float(da.comment[0]),4)
da.temperature = temperature
da.wavelength = wavelength
# cut the zero at the end
for i in range(len(data)):
data[i] = data[i,:,:-1]
return data
resoe = readEMUdat(js.examples.datapath + '/sample_5K_Q.dat.gz', temperature=5)
reso = js.dynamic.convert_e2w(resoe, 0, unit='meV')
# convert resolution to normalized data that elastic has meaning
for va in reso:
va.Y = va.Y/np.trapezoid(va.Y,va.X) # normalize resolution
# read new data and convert to 1/ns
same = readEMUdat(js.examples.datapath + '/sample_413K_Q.dat.gz',temperature=413)
sam = js.dynamic.convert_e2w(same, 293, unit='meV')
# Fit one after the other
for n in np.r_[0:len(sam)]:
#sam[n].makeErrPlot(yscale='log', title=str(sam[n].Q) + ' A\S-1')
sam[n].setlimit() # removes limits
sam[n].setlimit(elastic=[0, 1], bgr=[0.0001, 0.01]) # bgr lower bound arbitrary, should be >0
sam[n].setlimit(ampl=[0,10],amp2=[0]) # no negative amplitudes
sam[n].setlimit(g1=[0,100],g2=[0,100]) # no negative and some max fwhm
sam[n].setlimit(beta1=[0,1],beta2=[0,1]) # no negative and some max fwhm
sam[n].fit(js.dynamic.doubleStretchedExp_w,
{'elastic': 1, 'bgr': 2e-4, 'amp1': 0.2,'g1':16, 'beta1':1},
{'wmax':8, 'amp2':0,'g2': 1,'resolution':reso},
{'w': 'X'},
method='lm',max_nfev=20000)
#sam[n].killErrPlot() # close it , or uncomment to have them open
# look at one single errplot
n=3
sam[n].showlastErrPlot(title=str(sam[n].Q) , yscale='log')
sam[n].errPlot(sam[n].lastfit.X,sam[n].lastfit._se1,sy=0,li=[1,2,3],le='dexp 1') # add
sam[n].errPlot(sam[n].lastfit.X,sam[n].lastfit._se2,sy=0,li=[1,2,4],le='dexp 2')
sam[n].errPlot(sam[n].lastfit.X,sam[n].lastfit._el,sy=0,li=[1,2,5],le='elastic')
# sam[n].errplot.save(js.examples.imagepath+'/doubleStretchedExp_w.jpg', size=(2,2))
.. image:: ../../examples/images/doubleStretchedExp_w.jpg
:align: center
:width: 50 %
:alt: dynamic_t2f_examples.
"""
if wmax is None:
wmax = max(np.abs(w))
# get right Q value of resolution for selected vana and prune to max frequency
resolutionQ = resolution.filter(Q=Q)[0]
resolutionQwmax = resolutionQ.prune(-wmax,wmax,weight=None)
# get larger W range to remove spectral leakage
ddw = np.diff(w).mean()
nw = int(wmax/ddw)
ww = np.r_[w.min() - ddw * np.r_[nw:0:-1],w,w.max() + ddw * np.r_[1:nw+1]]
# two Lorentz model and elastic , all same w , Lorentz is normalized to area=1
model1 = stretchedExp_w(ww, g1, beta1)
model1.Y = amp1* model1.Y
# model1 = transDiff_w(ww, q=Q, D=fwhm1) #
model2 = stretchedExp_w(ww, g2, beta2)
model2.Y = amp2* model2.Y
# sum amplitudes
both = model1.copy()
both.Y = model1.Y + model2.Y
# here interpolate more cuts the extension of convolve that happens because of prune
convboth = convolve(both, resolutionQwmax, normB=True).interpolate(w)
convmodel1 = convolve(model1, resolutionQwmax, normB=True).interpolate(w)
convmodel2 = convolve(model2, resolutionQwmax, normB=True).interpolate(w)
# compose the data to return
# use resolution as elastic contribution
result = dA(np.c_[w, convboth.Y+elastic*resolutionQ.Y+bgr,
convmodel1.Y+bgr,
convmodel2.Y+bgr,
elastic*resolutionQ.Y+bgr].T)
# append some of the parameters
result.columnname='w;all;se1;se2;el' # name the columns for easier access
result.elastic = elastic
result.g1 = g1
result.g2 = g2
result.amp1 =amp1
result.amp2 =amp2
result.setColumnIndex(iey=None)
result.modelname = inspect.currentframe().f_code.co_name
return result
[docs]
def threeLorentz_w(w, Q, g1, a1, resolution, g2=1, a2=0, g3=1, a3=0, a0=0, bgr=0, u2 = 0, wmax=None):
r"""
Three Lorents functions with elastic contribution and resolution smearing.
A convenience function that allows to directly fit experimental data.
Parameters
----------
w : array
Frequencies 1/ns
Q : float
Wavevector in 1/nm for selecting the correct resolution.
g1,g2, g3 : float
FWHM units 1/ns.
a1,a2,a3 : float
Amplitudes of the Lorentz functions.
a0 : float
Amplitude elastic contribution. If the resolution is normalized it is directly the elastic contribution.
resolution : dataList
Measured or calculated resolution function with same Q as in the data to fit.
- Normalize the resolution before usage that the elastic amplitude a0 gets meaningful.
- In general it is assumed that the resolution and data w are equidistant.
- A fitted resolution can be used to smooth noise in a measured resolution.
See :py:func:`~.dynamic.frequencydomain.resolution_w` how to fit and use the '.lastfit' .
bgr : float
Background.
u2 : float , default=0
Means square displacement :math:`\langle u^2\rangle` for the Debye-Waller factor in 1/Q**2 units.
Only useful if several Q are fitted simultaneous.
The used Debye-Waller factor is :math:`DWF(Q,\langle u^2\rangle) = exp(-Q^2\langle u^2\rangle /3)`
wmax : float, default None
Cutoff frequency for resolution. Not needed for fitted resolution.
Returns
-------
model : dataArray
Q : wavevector
Notes
-----
We sum here three Lorentz functions with amplitudes and add an elastic contribution.
For Lorentz see :py:func:`lorentz_w` that are convoluted with the resolution cut at wmax.
As elastic contribution we use the given resolution cut at wmax and multiplied by value `elastic`.
For independent process j the inelastic measurement is often expressed by an elastic and an inelastic term as
.. math:: S_j(Q,\omega) \propto [B_{0j}(Q)\delta(Q) + (1-B_{0j}(Q)) L_j(Q,\omega)] \otimes R(Q,\omega)
with the resolution function :math:`R(Q, \omega)` and the Lorentzian
.. math:: L_j(Q,\omega) = \frac{1}{\pi} \frac{\Gamma_j(Q)}{\omega^2+\Gamma^2_j(Q)}
For several processes the summed intensity including a Debye-Waller factor is
.. math:: I(Q,\omega) = [ exp(-Q^2\langle u^2\rangle) \big( &A_0(Q)\delta(Q) + \\
\sum_j &A_j(Q) L_j(Q,\omega) \big) ] \otimes R(Q,\omega)
where the elastic contributions are summed in a single term :math:`A_{0}(Q)\delta(Q)`
Examples
--------
::
import jscatter as js
import numpy as np
def readEMUdat(filename, wavelength=6.27084, temperature=0):
# Read data form EMU@ANSTO as direct output from Mantid (No .inx converison needed)
data = js.dL(filename, delimiter=',')
for da in data:
da.Q = np.round(float(da.comment[0]),4)
da.temperature = temperature
da.wavelength = wavelength
# cut the zero at the end
for i in range(len(data)):
data[i] = data[i,:,:-1]
return data
vanae = readEMUdat(js.examples.datapath + '/sample_5K_Q.dat.gz', temperature=5)
vana = js.dynamic.convert_e2w(vanae, 0, unit='meV')
# convert resolution to normalized data that elastic has meaning
for va in vana:
va.integral = np.trapezoid(va.Y,va.X)
va.Y = va.Y/va.integral # normalise resolution
va.eY = va.eY/va.integral # normalise resolution
# read new data and convert to 1/ns
same = readEMUdat(js.examples.datapath + '/sample_413K_Q.dat.gz',temperature=413)
sam = js.dynamic.convert_e2w(same, 293, unit='meV')
# Fit one after the other
for n in np.r_[0:len(sam)]:
#sam[n].makeErrPlot(yscale='log', title=str(sam[1].Q) + ' A\S-1')
sam[n].setlimit() # removes limits
sam[n].setlimit(a0=[0, 1], bgr=[0., 0.01])
sam[n].setlimit(a1=[0],a2=[0],a3=[0]) # no negative amplitudes
sam[n].setlimit(g1=[0,100],g2=[0,100],g3=[0,100]) # no negative and some max
sam[n].fit(js.dynamic.threeLorentz_w,
{'a0': 1, 'bgr': 1e-5, 'a1': 0.2,'g1':16,'a2': 0.2,'g2':6, },
{'wmax':8, 'a3':0,'g3': 1,'resolution':vana},
{'w': 'X'},
method='lm',max_nfev=20000)
#sam[n].killErrPlot() # close it , or uncomment to have them open
# look at one single errplot
n=3
sam[n].showlastErrPlot(title=str(sam[n].Q) , yscale='log')
sam[n].errPlot(sam[n].lastfit.X,sam[n].lastfit._lo1,sy=0,li=[1,2,3],le='lo 1') # add
sam[n].errPlot(sam[n].lastfit.X,sam[n].lastfit._lo2,sy=0,li=[1,2,4],le='lo 2')
sam[n].errPlot(sam[n].lastfit.X,sam[n].lastfit._lo3,sy=0,li=[1,2,6],le='lo 3')
sam[n].errPlot(sam[n].lastfit.X,sam[n].lastfit._el,sy=0,li=[1,2,5],le='elastic')
# sam[n].errplot.save(js.examples.imagepath+'/threeLorentz_w.jpg', size=(2,2))
.. image:: ../../examples/images/threeLorentz_w.jpg
:align: center
:width: 50 %
:alt: dynamic_t2f_examples.
"""
if wmax is None:
wmax = max(np.abs(w))
# get right Q value of resolution for selected vana and prune to max frequency
resolutionQ = resolution.filter(Q=Q)[0]
resolutionQwmax = resolutionQ.prune(-wmax, wmax, weight=None)
# get larger W range to remove spectral leakage
ddw = np.diff(w).mean()
nw = int(wmax / ddw)
ww = np.r_[w.min() - ddw * np.r_[nw:0:-1], w, w.max() + ddw * np.r_[1:nw + 1]]
# two Lorentz model and elastic , all same w , Lorentz is normalized to area=1
model1 = lorentz_w(ww, g1)
model1.Y = a1 * model1.Y
model2 = lorentz_w(ww, g2)
model2.Y = a2 * model2.Y
model3 = lorentz_w(ww, g3)
model3.Y = a3 * model2.Y
# sum amplitudes
all = model1.copy()
all.Y = model1.Y + model2.Y + model3.Y
# here interpolate more cuts the extension of convolve that happens because of prune
convall = convolve(all, resolutionQwmax, normB=True).interpolate(w)
convmodel1 = convolve(model1, resolutionQwmax, normB=True).interpolate(w)
convmodel2 = convolve(model2, resolutionQwmax, normB=True).interpolate(w)
convmodel3 = convolve(model3, resolutionQwmax, normB=True).interpolate(w)
# Debye Waller factor
if u2 >0:
DWF = np.exp(-Q**2*u2/3)
else:
DWF = 1
# compose the data to return
# use resolution as elastic contribution
result = dA(np.c_[w, DWF * (convall.Y + a0 * resolutionQ.Y) + bgr,
DWF * convmodel1.Y + bgr,
DWF * convmodel2.Y + bgr,
DWF * convmodel3.Y + bgr,
DWF * a0 * resolutionQ.Y + bgr].T)
# append some of the parameters
result.columnname = 'w;all;lo1;lo2;lo3;el' # name the columns for easier access
result.a0 = a0
result.g1 = g1
result.g2 = g2
result.g3 = g3
result.a1 = a1
result.a2 = a2
result.a3 = a3
result.DWF = DWF
result.u2 = u2
result.setColumnIndex(iey=None)
result.modelname = inspect.currentframe().f_code.co_name
return result
[docs]
def fixedWindowScanArrhenius(Q, Temp, w, u2, A1, B1, Ea_1, tau0_1,tau_res, bgr=0,
A2=0, B2=0, Ea_2=None, tau0_2=None,
A3=0, B3=0, Ea_3=None, tau0_3=None):
r"""
Intensities of fixed window scans assuming Arrhenius behavior of up to 3 processes for
elastic and inelastic FWS. Shortcut fWSA
In case of thermally activated hindered motions the characteristic time :math:`\tau` of a process
is represented by an Arrhenius ansatz
:math:`\tau = \tau_{\infty}e^{E_a / RT}` [1]_.
If the elastic FWS has contribution of the inelastic FWS it needs to be included in the bgr.
Parameters
----------
Temp : array, float
Temperature in Kelvin.
Q : array, float
Scattering vector in units 1/nm.
w : float
Offset frequency :math:`\omega_{off}` in units 1/ns.
Offset energy :math:`\Delta E[\mu eV] /\hbar= \omega_{off}` with ``js.dynamic.hbar`` .
u2 : array len(Q), float
Means square displacements for the Debye-Waller factor in 1/Q**2 units.
- float : assumes u2*Temp => <u2> with given value as scaling with temperature.
- dataAray : temperature dependence with Temp as X and u2 as Y .
Values u2(Temp) will be linear interpolated from dataArray.
B1,B2,B3 : float
Intensity of the j-th process :math:`B_j` related to the number of atoms contributing to it.
A1,A2,A3 : float
EISF, elastic fraction of the j-th dynamic process contributing to elastic FWS.
(1-Aj) is the inelastic part.
Ea_1,Ea_2,Ea_3 : float
Activation energy :math:`E_{aj}` of the j-th dynamic process in units meV.
tau0_1,tau0_2,tau0_3 : float
The high-temperature relaxation time limit of the j-th dynamic process
:math:`τ_{0j}(Q)` in 1/ns.
tau_res : float
1/HWHM of instrumental resolution in ns.
Can be extracted from a resolution measurement at least approximately
by :py:func:`~jscatter.dynamic.fft.getHWHM` or by a Lorentz fit to the resolution.
bgr : array len(Q)
For w=0 the elastic background, but for w>0 the background contribution from the
wings of the resolution function.
Returns
-------
: dataArray
columnname 'Temp;IFWS;IFWS1;IFWS2;IFWS3;DWF;u2'
Notes
-----
We use here up to n=3 processes similar to [3]_ equ 2 + 3.
If only a single measurement (one :math:`w_{off}`) is fitted some parameters have to be fixed (Bj=1):
**elastic FWS** :math:`w_{off}=0`:
.. math:: I_{e}(Q,T) = DWF(Q) \big[ bgr_{e} +
\sum_j^n B_j(Q) \Bigl\{
A_{0j}(Q) + \frac{2}{\pi} (1-A_{0j}(Q)) arctan(\frac{\tau_j(Q)}{\tau_{res}})
\Bigr\}
\big]
**inelastic FWS** :math:`w_{off}\neq 0`:
.. math:: I_{i}(Q,T) = DWF(Q) \big[
bgr_{i} + \sum_j^n B_j(Q) \frac{1}{\pi}(1-A_{0j}(Q)) \frac{\tau_j(Q)}{1+\omega_{off}^2\tau_j^2(Q)}
\big]
with :math:`A_{0j}(Q)` being the EISF of the j-th dynamic process, :math:`\tau_{res}`
is the instrumental resolution and the Debye-Waller factor :math:`DWF(Q)=exp(-Q^2 \langle u^2 \rangle /3)`.
.. math:: \tau_j(Q) = \tau_{0j}(Q)exp(\frac{E_{aj}}{k_BT})
is the relaxation time of the j-th dynamic process.
:math:`\tau(Q)` is related to the HWHM Γ of a Lorentzian function commonly employed to model QENS broadening
:math:`\Gamma(Q)=1/\tau(Q)`.
From [3]_:
The equations above rely on the assumption that the described dynamic processes follow the Arrhenius law,
which is not always the case. They also assume independent motion, which otherwise would require a
convolution of two processes if they influence one another. The resolution is accounted for only in
terms of its width, and in the case of the IFWS can still contribute to the background and
influence the scaling parameter B in eq 3 (->inelastic FWS).
Therefore, initial information can be obtained from the FWS only, yet these assumptions
need to be taken into account.
**Connection to full inelastic scans** :
For each process j the full inelastic measurement can be expressed as
.. math:: S_j(Q,\omega) = A_{0j}(Q)\delta(Q) + (1-A_{0j}(Q)L(Q,\omega))
with the Lorentzian
.. math:: L_j(Q,\omega) = \frac{1}{\pi} \frac{\tau_j(Q)}{1+\omega^2\tau^2_j(Q)}
The intensity with the relative contributions :math:`B_j(Q)` of the processes is :
.. math:: I(Q,\omega) = [DWF(Q)\sum_j B_j(Q) S_j(Q,\omega)] \otimes R(Q,\omega)
where the elastic contribution is summed :math:`A_{0}(Q)\delta(Q) = \sum_j A_{0j}(Q)\delta(Q)`
(see :py:func:`threeLorentz_w`). The resolution function is :math:`R(Q, \omega)`.
Examples
--------
Example data showing 2 processes.
First only one ``woff`` then both together
We allow here 2 activation energies ``Ea_1`` dependent of ``woff`` because of e.g. temperature up and down scans
resulting in different average temperatures.
Parameters ``A1`` and ``A2`` are the same for same Q while Bi are independent for all measurements.
Other parameters are for all the same.
If you know the Q dependence of Ai and Bi of a process these can be included by encapsulating this model into
a small wrapper function.
Important is that the ``woff``, ``tau_res`` are stored in the data.
If separate ``Q`` values are stored in different files with Q in filename just read all and
recover the ``Q`` from filename like ``woff``.
The model can be fitted together with full scans like was done in [3]_ (example will follow later).
::
import jscatter as js
import numpy as np
fixedWindowScan = js.dynamic.fixedWindowScanArrhenius
# data with woff =0 and 3 meV; Q is found in data
fws = js.dL(js.examples.datapath +'/FWSsun_?microeV.dat.gz')
for data in fws:
# add woff from name in µeV convert to 1/ns
data.woff = round(float(data.name.split('_')[1].split('mic')[0])/ js.dynamic.hbar,3)
# tau_res from resolution of instrument, could be Q dependent
data.tau_res = 1.66 # resolution width in ns
data.I0 = data.Y[:5].mean() # start value for B
# only elastic contribution with two process and independent I0
fws0 = fws.filter(woff=0)
fws0.makeErrPlot(yscale='log',residuals='relative')
# resonable limits
fws0.setlimit(u2=[0,1,0,2],bgr=[0],I0=[0,10])
fws0.setlimit(A1=[0,1], B1=[0,10],Ea_1=[50,500],tau0_1=[0,0.01])
fws0.setlimit(A2=[0,1], B2=[0,10],Ea_2=[0,50],tau0_2=[0,0.02])
fws0.setlimit(A3=[0,1], B3=[0,10],Ea_3=[50,500],tau0_3=[0,0.02])
# create reasonable starting values, use tolist to get list
I0 = fws0.I0.array
I0[I0<1] *=50
fws0.fit(model=fixedWindowScan,
freepar={'bgr':[0.004],'u2':0.0002,
'A1':[.15,'Q'], 'B1':I0.tolist(), 'Ea_1':270, 'tau0_1':2e-5,
'A2':[0.05,'Q'], 'B2':[0.03], 'Ea_2':30, 'tau0_2':1e-2,
},
fixpar={'A3':.05, 'B3':0, 'Ea_3':200, 'tau0_3':0.00001},
mapNames={'Temp':'X','w':'woff'},workers=0) # ,method='Nelder-Mead')
# EFWS + IFWS fit together with two process and independent I0
fws.makeErrPlot(yscale='log',residuals='relative')
# resonable limits
fws.setlimit(u2=[0,1,0,2],bgr=[0],I0=[0,10])
fws.setlimit(A1=[0,1], B1=[0,10],Ea_1=[50,500],tau0_1=[0,0.01])
fws.setlimit(A2=[0,1], B2=[0,10],Ea_2=[0,50],tau0_2=[0,0.02])
fws.setlimit(A3=[0,1], B3=[0,10],Ea_3=[50,500],tau0_3=[0,0.02])
# create reasonable starting values, use tolist to get list
I0 = fws.I0.array
I0[I0<1] *=50
fws.fit(model=fixedWindowScan,
freepar={'bgr':[0.004,0.001,'woff'],'u2':0.0002, #'I0':fws.I0,
'A1':[.15,'Q'], 'B1':I0.tolist(), 'Ea_1':[270,210,'woff'], 'tau0_1':2e-5,
'A2':[0.05,'Q'], 'B2':[0.03], 'Ea_2':30, 'tau0_2':1e-2,
},
fixpar={'A3':.05, 'B3':0, 'Ea_3':200, 'tau0_3':0.00001},
mapNames={'Temp':'X','w':'woff'},workers=0) # ,method='Nelder-Mead')
fws.errplot.legend(x=5,y=3,charsize=0.5)
# fws.savelastErrPlot(js.examples.imagepath+'/fixedWindoScanfit.jpg', size=(2,1),dpi=600)
.. image:: ../../examples/images/fixedWindoScanfit.jpg
:align: center
:width: 50 %
:alt: fixedWindoScanfit.
References
----------
.. [1] An investigation of micro-brownian motions in polydimethylsiloxane by complementary
incoherent-neutron-scattering
and nuclear-magnetic-resonance experiments below room temperature.
Grapengeter, H. H., Alefeld, B. & Kosfeld, R.
Colloid & Polymer Science 265, (1987). https://dx.doi.org/10.1007/BF01412711
.. [2] New Possibilities with Inelastic Fixed Window Scans and Linear Motor Doppler Drives on
High Resolution Neutron Backscattering Spectrometers.
Frick, B.; Combet, J.; Van Eijck, L.
Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and
Associated Equipment 2012, 669, 7–13. https://doi.org/10.1016/J.NIMA.2011.11.090.
.. [3] Quasi-Elastic Neutron Scattering of Citrate-Capped Iron Oxide Nanoparticles: Distinguishing between Ligand,
Water, and Magnetic Dynamics.
Plekhanov, M. S.; Thomä, S. L. J.; Magerl, A.; Appel, M.; Zobel, M.
J. Phys. Chem. C 2024, 128 (28), 11661–11671. https://doi.org/10.1021/acs.jpcc.4c00479.
"""
Q = np.atleast_1d(Q)
Temp = np.atleast_1d(Temp)
w = np.atleast_1d(w)
# determine u2 for DWF
if isinstance(u2, dA):
u2 = u2.interp(Temp)
else:
u2 = np.atleast_1d(u2) * Temp
# Debye Waller factor
DWF = np.exp(-Q ** 2 * u2)
IFWS2 = np.zeros_like(Temp)
IFWS3 = np.zeros_like(Temp)
kT = k * Temp
if w == 0:
# elastic fixed Window scan
tau_1 = tau0_1 * np.exp(np.clip(Ea_1 / kT, 0, 700))
IFWS1 = A1 + 2 / pi * (1 - A1) * np.arctan(tau_1 / tau_res)
if B2 > 0:
tau_2 = tau0_2 * np.exp(np.clip(Ea_2 / kT, 0, 700))
IFWS2 = A2 + 2 / pi * (1 - A2) * np.arctan(tau_2 / tau_res)
if B3 > 0:
tau_3 = tau0_3 * np.exp(np.clip(Ea_3 / kT, 0, 700))
IFWS3 = A3 + 2 / pi * (1 - A3) * np.arctan(tau_3 / tau_res)
result = dA(np.c_[Temp, DWF * (B1 * IFWS1 + B2 * IFWS2 + B3 * IFWS3 + bgr),
DWF * B1 * IFWS1,
DWF * B2 * IFWS2,
DWF * B3 * IFWS3,
DWF, u2].T)
result.type = 'elasticFWS'
else:
# inelastic fixed Window scan at offset w
# using inverse tau avoids quadratic overflow
taui = 1 / (tau0_1 * np.exp(np.clip(Ea_1 / kT, 0, 700)))
IFWS1 = (1 - A1) / pi * (taui / (taui ** 2 + w ** 2))
if B2 > 0:
taui = 1 / (tau0_2 * np.exp(np.clip(Ea_2 / kT, 0, 700)))
IFWS2 = (1 - A2) / pi * (taui / (taui ** 2 + w ** 2))
if B3 > 0:
taui = 1 / (tau0_3 * np.exp(np.clip(Ea_3 / kT, 0, 700)))
IFWS3 = (1 - A3) / pi * (taui / (taui ** 2 + w ** 2))
result = dA(np.c_[Temp, DWF * (B1 * IFWS1 + B2 * IFWS2 + B3 * IFWS3 + bgr),
DWF * B1 * IFWS1,
DWF * B2 * IFWS2,
DWF * B3 * IFWS3,
DWF, u2].T)
result.type = 'inelasticFWS'
result.columnname = 'Temp;FWS;FWS1;FWS2;FWS3;DWF;u2'
result.setColumnIndex(iey=None) # no errors in models
result.modelname = inspect.currentframe().f_code.co_name
result.w = w
result.energy = w
result.Q = Q
result.A1 = A1
result.B1 = B1
result.Ea_1 = Ea_1
result.tau0_1 = tau0_1
if B2 > 0:
result.A2 = A2
result.B2 = B2
result.Ea_2 = Ea_2
result.tau0_2 = tau0_2
if B3 > 0:
result.A3 = A3
result.B3 = B3
result.Ea_3 = Ea_3
result.tau0_3 = tau0_3
result.tau_res = tau_res
result.bgr = bgr
return result
fWSA = fixedWindowScanArrhenius