summaryrefslogtreecommitdiff
path: root/lmfit/lineshapes.py
diff options
context:
space:
mode:
Diffstat (limited to 'lmfit/lineshapes.py')
-rw-r--r--lmfit/lineshapes.py89
1 files changed, 46 insertions, 43 deletions
diff --git a/lmfit/lineshapes.py b/lmfit/lineshapes.py
index 064a18c..7ecdc42 100644
--- a/lmfit/lineshapes.py
+++ b/lmfit/lineshapes.py
@@ -1,25 +1,24 @@
-"""Basic model line shapes and distribution functions."""
+"""Basic model lineshapes and distribution functions."""
-import warnings
-
-from numpy import (arctan, copysign, cos, exp, isnan, log, pi, real, sin, sqrt,
- where)
+from numpy import (arctan, copysign, cos, exp, isclose, isnan, log, log1p,
+ maximum, minimum, pi, real, sin, sqrt, where)
+from scipy.special import betaln as betalnfcn
from scipy.special import erf, erfc
from scipy.special import gamma as gamfcn
+from scipy.special import loggamma as loggammafcn
from scipy.special import wofz
log2 = log(2)
s2pi = sqrt(2*pi)
-spi = sqrt(pi)
s2 = sqrt(2.0)
# tiny had been numpy.finfo(numpy.float64).eps ~=2.2e16.
# here, we explicitly set it to 1.e-15 == numpy.finfo(numpy.float64).resolution
tiny = 1.0e-15
functions = ('gaussian', 'gaussian2d', 'lorentzian', 'voigt', 'pvoigt',
- 'moffat', 'pearson7', 'breit_wigner', 'damped_oscillator',
+ 'moffat', 'pearson4', 'pearson7', 'breit_wigner', 'damped_oscillator',
'dho', 'logistic', 'lognormal', 'students_t', 'expgaussian',
- 'doniach', 'donaich', 'skewed_gaussian', 'skewed_voigt',
+ 'doniach', 'skewed_gaussian', 'skewed_voigt',
'thermal_distribution', 'step', 'rectangle', 'exponential',
'powerlaw', 'linear', 'parabolic', 'sine', 'expsine',
'split_lorentzian')
@@ -148,6 +147,26 @@ def moffat(x, amplitude=1, center=0., sigma=1, beta=1.):
return amplitude / (((x - center)/max(tiny, sigma))**2 + 1)**beta
+def pearson4(x, amplitude=1.0, center=0.0, sigma=1.0, expon=1.0, skew=0.0):
+ """Return a Pearson4 lineshape.
+
+ Using the Wikipedia definition:
+
+ pearson4(x, amplitude, center, sigma, expon, skew) =
+ amplitude*|gamma(expon + I skew/2)/gamma(m)|**2/(w*beta(expon-0.5, 0.5)) * (1+arg**2)**(-expon) * exp(-skew * arctan(arg))
+
+ where ``arg = (x-center)/sigma``, `gamma` is the gamma function and `beta` is the beta function.
+
+ For more information, see: https://en.wikipedia.org/wiki/Pearson_distribution#The_Pearson_type_IV_distribution
+
+ """
+ expon = max(tiny, expon)
+ sigma = max(tiny, sigma)
+ arg = (x - center) / sigma
+ logprefactor = 2 * (real(loggammafcn(expon + skew * 0.5j)) - loggammafcn(expon)) - betalnfcn(expon - 0.5, 0.5)
+ return (amplitude / sigma) * exp(logprefactor - expon * log1p(arg * arg) - skew * arctan(arg))
+
+
def pearson7(x, amplitude=1.0, center=0.0, sigma=1.0, expon=1.0):
"""Return a Pearson7 lineshape.
@@ -188,7 +207,7 @@ def damped_oscillator(x, amplitude=1.0, center=1., sigma=0.1):
return amplitude/sqrt((1.0 - (x/center)**2)**2 + (2*sigma*x/center)**2)
-def dho(x, amplitude=1., center=0., sigma=1., gamma=1.0):
+def dho(x, amplitude=1., center=1., sigma=1., gamma=1.0):
"""Return a Damped Harmonic Oscillator.
Similar to the version from PAN:
@@ -201,16 +220,19 @@ def dho(x, amplitude=1., center=0., sigma=1., gamma=1.0):
``lp(x, center, sigma) = 1.0 / ((x+center)**2 + sigma**2)``
"""
+ factor = amplitude * sigma / pi
bose = (1.0 - exp(-x/max(tiny, gamma)))
if isinstance(bose, (int, float)):
- bose = max(tiny, bose)
+ bose = not_zero(bose)
else:
bose[where(isnan(bose))] = tiny
- bose[where(bose <= tiny)] = tiny
+ bose[where(abs(bose) <= tiny)] = tiny
lm = 1.0/((x-center)**2 + sigma**2)
lp = 1.0/((x+center)**2 + sigma**2)
- return amplitude*sigma/pi*(lm - lp)/bose
+ return factor * where(isclose(x, 0.0),
+ 4*gamma*center/(center**2+sigma**2)**2,
+ (lm - lp)/bose)
def logistic(x, amplitude=1., center=0., sigma=1.):
@@ -289,20 +311,6 @@ def doniach(x, amplitude=1.0, center=0, sigma=1.0, gamma=0.0):
return scale*cos(pi*gamma/2 + gm1*arctan(arg))/(1 + arg**2)**(gm1/2)
-def donaich(x, amplitude=1.0, center=0, sigma=1.0, gamma=0.0):
- """Return a Doniach Sunjic asymmetric lineshape.
-
- Function added here for backwards-compatibility, will emit a
- `FutureWarning` when used.
-
- """
- msg = ('Please correct the name of your lineshape function: donaich --> '
- 'doniach. The incorrect spelling will be removed in a later '
- 'release.')
- warnings.warn(FutureWarning(msg))
- return doniach(x, amplitude, center, sigma, gamma)
-
-
def skewed_gaussian(x, amplitude=1.0, center=0.0, sigma=1.0, gamma=0.0):
"""Return a Gaussian lineshape, skewed with error function.
@@ -397,8 +405,8 @@ def thermal_distribution(x, amplitude=1.0, center=0.0, kt=1.0, form='bose'):
elif form.startswith('fermi'):
offset = 1
else:
- msg = "Invalid value ('%s') for argument 'form'; should be one of %s."\
- % (form, "'maxwell', 'fermi', or 'bose'")
+ msg = (f"Invalid value ('{form}') for argument 'form'; should be one "
+ "of 'maxwell', 'fermi', or 'bose'.")
raise ValueError(msg)
return real(1/(amplitude*exp((x - center)/not_zero(kt)) + offset + tiny*1j))
@@ -410,7 +418,7 @@ def step(x, amplitude=1.0, center=0.0, sigma=1.0, form='linear'):
Starts at 0.0, ends at `amplitude`, with half-max at `center`, and
rising with `form`:
- - `'linear'` (default) = amplitude * min(1, max(0, arg))
+ - `'linear'` (default) = amplitude * min(1, max(0, arg + 0.5))
- `'atan'`, `'arctan'` = amplitude * (0.5 + atan(arg)/pi)
- `'erf'` = amplitude * (1 + erf(arg))/2.0
- `'logistic'` = amplitude * [1 - 1/(1 + exp(arg))]
@@ -423,15 +431,14 @@ def step(x, amplitude=1.0, center=0.0, sigma=1.0, form='linear'):
if form == 'erf':
out = 0.5*(1 + erf(out))
elif form == 'logistic':
- out = (1. - 1./(1. + exp(out)))
+ out = 1. - 1./(1. + exp(out))
elif form in ('atan', 'arctan'):
out = 0.5 + arctan(out)/pi
elif form == 'linear':
- out[where(out < 0)] = 0.0
- out[where(out > 1)] = 1.0
+ out = minimum(1, maximum(0, out + 0.5))
else:
- msg = "Invalid value ('%s') for argument 'form'; should be one of %s."\
- % (form, "'erf', 'logistic', 'atan', 'arctan', or 'linear'")
+ msg = (f"Invalid value ('{form}') for argument 'form'; should be one "
+ "of 'erf', 'logistic', 'atan', 'arctan', or 'linear'.")
raise ValueError(msg)
return amplitude*out
@@ -462,18 +469,14 @@ def rectangle(x, amplitude=1.0, center1=0.0, sigma1=1.0,
if form == 'erf':
out = 0.5*(erf(arg1) + erf(arg2))
elif form == 'logistic':
- out = (1. - 1./(1. + exp(arg1)) - 1./(1. + exp(arg2)))
+ out = 1. - 1./(1. + exp(arg1)) - 1./(1. + exp(arg2))
elif form in ('atan', 'arctan'):
out = (arctan(arg1) + arctan(arg2))/pi
elif form == 'linear':
- arg1[where(arg1 < 0)] = 0.0
- arg1[where(arg1 > 1)] = 1.0
- arg2[where(arg2 > 0)] = 0.0
- arg2[where(arg2 < -1)] = -1.0
- out = arg1 + arg2
+ out = 0.5*(minimum(1, maximum(-1, arg1)) + minimum(1, maximum(-1, arg2)))
else:
- msg = "Invalid value ('%s') for argument 'form'; should be one of %s."\
- % (form, "'erf', 'logistic', 'atan', 'arctan', or 'linear'")
+ msg = (f"Invalid value ('{form}') for argument 'form'; should be one "
+ "of 'erf', 'logistic', 'atan', 'arctan', or 'linear'.")
raise ValueError(msg)
return amplitude*out
@@ -501,7 +504,7 @@ def powerlaw(x, amplitude=1, exponent=1.0):
def linear(x, slope=1.0, intercept=0.0):
"""Return a linear function.
- linear(x, slope, interceps) = slope * x + intercept
+ linear(x, slope, intercept) = slope * x + intercept
"""
return slope * x + intercept
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-26 01:48:57 +0000