Components¶

Components are the basic ingredients for all models. Currently, the available components are Gaussian, Lorentzian, Voigt (the convolution of a Gaussian with a Lorentzian), delta function, damped harmonic oscillator and polynomial. This notebooks shows how to use the components.

Note in particular that a Gaussian, Lorentzian, Voigt or delta function where the center has not been given will be centered at 0.

In [1]:

import matplotlib.pyplot as plt
import numpy as np
import scipp as sc

import easydynamics.sample_model as sm

%matplotlib widget

import matplotlib.pyplot as plt import numpy as np import scipp as sc import easydynamics.sample_model as sm %matplotlib widget

In [2]:

# Creating a component
gaussian = sm.Gaussian(display_name='Gaussian', width=0.5, area=1)
dho = sm.DampedHarmonicOscillator(display_name='DHO', center=1.0, width=0.3, area=2.0)
lorentzian = sm.Lorentzian(display_name='Lorentzian', center=-1.0, width=0.2, area=1.0)
polynomial = sm.Polynomial(
    display_name='Polynomial', coefficients=[-0.2, 0, 0.5]
)  # y=-0.2+0.5*x^2
exponential = sm.Exponential(display_name='Exponential', amplitude=1.0, rate=-0.5)

x = np.linspace(-2, 2, 100)

plt.figure()
y = gaussian.evaluate(x)
plt.plot(x, y, label='Gaussian')
y = dho.evaluate(x)
plt.plot(x, y, label='DHO')
y = lorentzian.evaluate(x)
plt.plot(x, y, label='Lorentzian')
y = polynomial.evaluate(x)
plt.plot(x, y, label='Polynomial')
y = exponential.evaluate(x)
plt.plot(x, y, label='Exponential')
plt.legend()
plt.show()

# Creating a component gaussian = sm.Gaussian(display_name='Gaussian', width=0.5, area=1) dho = sm.DampedHarmonicOscillator(display_name='DHO', center=1.0, width=0.3, area=2.0) lorentzian = sm.Lorentzian(display_name='Lorentzian', center=-1.0, width=0.2, area=1.0) polynomial = sm.Polynomial( display_name='Polynomial', coefficients=[-0.2, 0, 0.5] ) # y=-0.2+0.5*x^2 exponential = sm.Exponential(display_name='Exponential', amplitude=1.0, rate=-0.5) x = np.linspace(-2, 2, 100) plt.figure() y = gaussian.evaluate(x) plt.plot(x, y, label='Gaussian') y = dho.evaluate(x) plt.plot(x, y, label='DHO') y = lorentzian.evaluate(x) plt.plot(x, y, label='Lorentzian') y = polynomial.evaluate(x) plt.plot(x, y, label='Polynomial') y = exponential.evaluate(x) plt.plot(x, y, label='Exponential') plt.legend() plt.show()

/tmp/ipykernel_2677/2610271038.py:19: UserWarning: The Polynomial with unique_name Polynomial_0 has negative values, which may not be physically meaningful.
  y = polynomial.evaluate(x)

In [3]:

# Suppress the warning from the Polynomial:
polynomial.suppress_warnings = True
y = polynomial.evaluate(x)

# Suppress the warning from the Polynomial: polynomial.suppress_warnings = True y = polynomial.evaluate(x)

In [4]:

# The area under the DHO curve is indeed equal to the area parameter.
xx = np.linspace(-15, 15, 10000)
yy = dho.evaluate(xx)
area = np.trapezoid(yy, xx)
print(f'Area under DHO curve: {area:.4f}')

# The area under the DHO curve is indeed equal to the area parameter. xx = np.linspace(-15, 15, 10000) yy = dho.evaluate(xx) area = np.trapezoid(yy, xx) print(f'Area under DHO curve: {area:.4f}')

Area under DHO curve: 1.9999

In [5]:

delta = sm.DeltaFunction(display_name='Delta', center=0.0, area=1.0)
x1 = np.linspace(-2, 2, 100)
y = delta.evaluate(x1)
x2 = np.linspace(-2, 2, 51)
y2 = delta.evaluate(x2)
plt.figure()
plt.plot(x1, y, label='Delta Function')
plt.plot(x2, y2, label='Delta Function (coarser)')
plt.legend()
plt.show()
# The area under the Delta function is indeed equal to
# the area parameter.
xx = np.linspace(-2, 2, 10000)
yy = delta.evaluate(xx)
area = np.trapezoid(y, x1)
print(area)

delta = sm.DeltaFunction(display_name='Delta', center=0.0, area=1.0) x1 = np.linspace(-2, 2, 100) y = delta.evaluate(x1) x2 = np.linspace(-2, 2, 51) y2 = delta.evaluate(x2) plt.figure() plt.plot(x1, y, label='Delta Function') plt.plot(x2, y2, label='Delta Function (coarser)') plt.legend() plt.show() # The area under the Delta function is indeed equal to # the area parameter. xx = np.linspace(-2, 2, 10000) yy = delta.evaluate(xx) area = np.trapezoid(y, x1) print(area)

0.9999999999999999

In [6]:

x1 = sc.linspace(dim='x', start=-2.0, stop=2.0, num=100, unit='meV')
x2 = sc.linspace(dim='x', start=-2.0 * 1e3, stop=2.0 * 1e3, num=101, unit='microeV')

polynomial = sm.Polynomial(display_name='Polynomial', coefficients=[0.1, 0, 0.5])  # y=0.1+0.5*x^2
y1 = polynomial.evaluate(x1)
y2 = polynomial.evaluate(x2)

plt.figure()
plt.plot(x1.values, y1, label='Polynomial meV', color='blue')
plt.plot(x2.values / 1000, y2, label='Polynomial microeV', linestyle='dashed', color='orange')
plt.legend()
plt.show()

x1 = sc.linspace(dim='x', start=-2.0, stop=2.0, num=100, unit='meV') x2 = sc.linspace(dim='x', start=-2.0 * 1e3, stop=2.0 * 1e3, num=101, unit='microeV') polynomial = sm.Polynomial(display_name='Polynomial', coefficients=[0.1, 0, 0.5]) # y=0.1+0.5*x^2 y1 = polynomial.evaluate(x1) y2 = polynomial.evaluate(x2) plt.figure() plt.plot(x1.values, y1, label='Polynomial meV', color='blue') plt.plot(x2.values / 1000, y2, label='Polynomial microeV', linestyle='dashed', color='orange') plt.legend() plt.show()

The ExpressionComponent can create almost any function you would like. Units are supported and must be given for each parameter either at construction as below, or afterwards as properties. The individual parameters can be accessed by their name in the usual way. You must use x fore the x-axis. It will warn if the units are not consistent.

In [7]:

expr = sm.ExpressionComponent(
    'A * exp(-(x - x0)**2 / (2*sigma**2)) +B*sin(2*pi*x/period)',
    parameters={'A': 10, 'x0': 0, 'sigma': 1},
    parameter_units={
        'A': 'dimensionless',
        'x0': 'meV',
        'sigma': 'meV',
        'B': 'dimensionless',
        'period': 'meV',
    },
    x_unit='meV',
    y_unit='dimensionless',
)

expr.A = 5
expr.sigma = 0.5

expr.period = 2.0

x = np.linspace(-5, 5, 100)
y = expr.evaluate(x)

plt.figure()
plt.plot(x, y, label='Expression Component')
plt.legend()
plt.show()
plt.xlabel('Energy (meV)')
plt.ylabel('Intensity (arb. units)')

expr = sm.ExpressionComponent( 'A * exp(-(x - x0)**2 / (2*sigma**2)) +B*sin(2*pi*x/period)', parameters={'A': 10, 'x0': 0, 'sigma': 1}, parameter_units={ 'A': 'dimensionless', 'x0': 'meV', 'sigma': 'meV', 'B': 'dimensionless', 'period': 'meV', }, x_unit='meV', y_unit='dimensionless', ) expr.A = 5 expr.sigma = 0.5 expr.period = 2.0 x = np.linspace(-5, 5, 100) y = expr.evaluate(x) plt.figure() plt.plot(x, y, label='Expression Component') plt.legend() plt.show() plt.xlabel('Energy (meV)') plt.ylabel('Intensity (arb. units)')

Out[7]:

Text(43.722222222222214, 0.5, 'Intensity (arb. units)')

In [8]:

expr = sm.ExpressionComponent(
    'A*erf(B*x)',
)

expr.A = 1.0
expr.B = 0.5


x = np.linspace(-5, 5, 100)
y = expr.evaluate(x)

plt.figure()
plt.plot(x, y, label='erf')
plt.legend()
plt.show()

expr = sm.ExpressionComponent( 'A*erf(B*x)', ) expr.A = 1.0 expr.B = 0.5 x = np.linspace(-5, 5, 100) y = expr.evaluate(x) plt.figure() plt.plot(x, y, label='erf') plt.legend() plt.show()