Skip to content

Kinetic model functions

Model functions are vectorised callables: f(T, Im, E, Tm, ...) -> ndarray. Obtain any model at runtime with tl.get_model("fo_rq").

First-order (FO)

tldecpy.models.fo.fo_rq

fo_rq(T, Im, E, Tm)

Evaluate first-order model fo_rq using Bos polynomial-rational :math:xQ(x).

Parameters:

Name Type Description Default
T ArrayLike

Temperature values in kelvin.

required
Im float

Peak intensity at :math:T_m.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature in kelvin.

required

Returns:

Type Description
ndarray

Intensity values :math:I(T).

References

.. [1] Bos, A. J. J., et al. (1993). Intercomparison synthetic glow curves (GLOCANIN). .. [2] Kitis, G., Gomes-Ros, J. M., and Tuyn, J. W. N. (1998). Thermoluminescence glow-curve deconvolution functions for first, second and general orders of kinetics.

Source code in tldecpy/models/fo.py
def fo_rq(T: ArrayLike, Im: float, E: float, Tm: float) -> FloatArray:
    r"""
    Evaluate first-order model ``fo_rq`` using Bos polynomial-rational :math:`xQ(x)`.

    Parameters
    ----------
    T : ArrayLike
        Temperature values in kelvin.
    Im : float
        Peak intensity at :math:`T_m`.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature in kelvin.

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.

    References
    ----------
    .. [1] Bos, A. J. J., et al. (1993). Intercomparison synthetic glow curves (GLOCANIN).
    .. [2] Kitis, G., Gomes-Ros, J. M., and Tuyn, J. W. N. (1998).
           Thermoluminescence glow-curve deconvolution functions for first,
           second and general orders of kinetics.
    """
    return _fo_core(
        T=T,
        Im=Im,
        E=E,
        Tm=Tm,
        Q_fn=lambda x: np.asarray(1.0 - alpha_poly(x), dtype=np.float64),
    )

tldecpy.models.fo.fo_rb

fo_rb(T, Im, E, Tm)

Evaluate first-order model fo_rb using AAA approximation for :math:Q(x).

Parameters:

Name Type Description Default
T ArrayLike

Temperature values in kelvin.

required
Im float

Peak intensity at :math:T_m.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature in kelvin.

required

Returns:

Type Description
ndarray

Intensity values :math:I(T).

Source code in tldecpy/models/fo.py
def fo_rb(T: ArrayLike, Im: float, E: float, Tm: float) -> FloatArray:
    r"""
    Evaluate first-order model ``fo_rb`` using AAA approximation for :math:`Q(x)`.

    Parameters
    ----------
    T : ArrayLike
        Temperature values in kelvin.
    Im : float
        Peak intensity at :math:`T_m`.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature in kelvin.

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.
    """
    try:
        return _fo_core(T=T, Im=Im, E=E, Tm=Tm, Q_fn=Q_aaa)
    except (FileNotFoundError, ModuleNotFoundError):
        warnings.warn(
            "AAA constants not found for fo_rb; falling back to (1 - alpha_poly) backend.",
            RuntimeWarning,
            stacklevel=2,
        )
        return _fo_core(
            T=T,
            Im=Im,
            E=E,
            Tm=Tm,
            Q_fn=lambda x: np.asarray(1.0 - alpha_poly(x), dtype=np.float64),
        )

tldecpy.models.fo.fo_ka

fo_ka(T, Im, E, Tm)

Evaluate first-order variant fo_ka (approximate closed form).

Parameters:

Name Type Description Default
T ArrayLike

Temperature values in kelvin.

required
Im float

Peak intensity at :math:T_m.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature in kelvin.

required

Returns:

Type Description
ndarray

Intensity values :math:I(T).

Source code in tldecpy/models/fo.py
def fo_ka(T: ArrayLike, Im: float, E: float, Tm: float) -> FloatArray:
    r"""
    Evaluate first-order variant ``fo_ka`` (approximate closed form).

    Parameters
    ----------
    T : ArrayLike
        Temperature values in kelvin.
    Im : float
        Peak intensity at :math:`T_m`.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature in kelvin.

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.
    """
    t_arr = np.maximum(np.asarray(T, dtype=np.float64), 1e-5)

    delta = (2.0 * KB_EV * t_arr) / E
    delta_m = (2.0 * KB_EV * Tm) / E
    arg = (E / (KB_EV * t_arr)) * ((t_arr - Tm) / Tm)
    exp_arg = np.exp(np.clip(arg, -50.0, 50.0))
    term = (t_arr**2 / Tm**2) * exp_arg * (1.0 - delta)
    exponent = 1.0 + arg - term - delta_m
    return np.asarray(Im * np.exp(np.clip(exponent, -100.0, 100.0)), dtype=np.float64)

tldecpy.models.fo.fo_wp

fo_wp(T, Im, E, Tm)

Evaluate first-order variant fo_wp (Weibull approximation).

Parameters:

Name Type Description Default
T ArrayLike

Temperature values in kelvin.

required
Im float

Peak intensity at :math:T_m.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature in kelvin.

required

Returns:

Type Description
ndarray

Intensity values :math:I(T).

Source code in tldecpy/models/fo.py
def fo_wp(T: ArrayLike, Im: float, E: float, Tm: float) -> FloatArray:
    r"""
    Evaluate first-order variant ``fo_wp`` (Weibull approximation).

    Parameters
    ----------
    T : ArrayLike
        Temperature values in kelvin.
    Im : float
        Peak intensity at :math:`T_m`.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature in kelvin.

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.
    """
    t_arr = np.asarray(T, dtype=np.float64)
    omega_approx = (2.52 * KB_EV * Tm**2) / (E + 2.0 * KB_EV * Tm)
    b_weibull = 6.44 * omega_approx
    x_val = (t_arr - Tm) / b_weibull + 0.996

    mask = x_val > 0.0
    y = np.zeros_like(t_arr, dtype=np.float64)
    x_safe = x_val[mask]
    term1 = x_safe**15
    term2 = np.exp(np.clip(-(x_safe**16), -100.0, 100.0))
    y[mask] = 2.713 * Im * term1 * term2
    return y

Second-order (SO)

tldecpy.models.so.so_ks

so_ks(T, Im, E, Tm)

Evaluate the Kitis-standard second-order expression.

Parameters:

Name Type Description Default
T ndarray

Temperature grid :math:T in kelvin.

required
Im float

Peak intensity :math:I_m.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature :math:T_m in kelvin.

required

Returns:

Type Description
ndarray

Intensity values :math:I(T).

References

.. [1] Kitis, G., Gomes-Ros, J. M., and Tuyn, J. W. N. (1998). Thermoluminescence glow-curve deconvolution functions for first, second and general orders of kinetics.

Source code in tldecpy/models/so.py
def so_ks(T: np.ndarray, Im: float, E: float, Tm: float) -> np.ndarray:
    r"""
    Evaluate the Kitis-standard second-order expression.

    Parameters
    ----------
    T : numpy.ndarray
        Temperature grid :math:`T` in kelvin.
    Im : float
        Peak intensity :math:`I_m`.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature :math:`T_m` in kelvin.

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.

    References
    ----------
    .. [1] Kitis, G., Gomes-Ros, J. M., and Tuyn, J. W. N. (1998).
           Thermoluminescence glow-curve deconvolution functions for first,
           second and general orders of kinetics.
    """
    temperature = np.maximum(np.asarray(T, dtype=float), 1e-5)
    k = KB_EV

    delta = (2.0 * k * temperature) / E
    delta_m = (2.0 * k * Tm) / E

    arg = (E / (k * temperature)) * ((temperature - Tm) / Tm)
    exp_arg = np.exp(np.clip(arg, -50.0, 50.0))

    numerator = 4.0 * Im * exp_arg
    denominator_term = (temperature**2 / Tm**2) * (1.0 - delta) * exp_arg + 1.0 + delta_m

    return numerator * (denominator_term**-2)

tldecpy.models.so.so_la

so_la(T, Im, E, Tm)

Evaluate the logistic-asymmetric second-order approximation.

Parameters:

Name Type Description Default
T ndarray

Temperature grid :math:T in kelvin.

required
Im float

Peak intensity :math:I_m.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature :math:T_m in kelvin.

required

Returns:

Type Description
ndarray

Intensity values :math:I(T).

References

.. [1] Pagonis, V., and Kitis, G. (2001). Fit of second-order thermoluminescence glow peaks using the logistic distribution function.

Source code in tldecpy/models/so.py
def so_la(T: np.ndarray, Im: float, E: float, Tm: float) -> np.ndarray:
    r"""
    Evaluate the logistic-asymmetric second-order approximation.

    Parameters
    ----------
    T : numpy.ndarray
        Temperature grid :math:`T` in kelvin.
    Im : float
        Peak intensity :math:`I_m`.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature :math:`T_m` in kelvin.

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.

    References
    ----------
    .. [1] Pagonis, V., and Kitis, G. (2001). Fit of second-order
           thermoluminescence glow peaks using the logistic distribution function.
    """
    temperature = np.asarray(T, dtype=float)
    k = KB_EV

    numerator = 1.189 * (Tm**4) * (k**2)
    denominator = E**2 + 4.0 * E * Tm * k
    a2 = np.sqrt(numerator / denominator)

    x_value = (temperature - Tm) / a2 + 0.38542
    exp_neg_x = np.exp(np.clip(-x_value, -50.0, 50.0))
    term = (1.0 + exp_neg_x) ** (-2.4702)
    return 5.2973 * Im * term * exp_neg_x

General-order (GO)

tldecpy.models.go.go_kg

go_kg(T, Im, E, Tm, b)

Evaluate the Kitis-general expression for general-order TL kinetics.

Parameters:

Name Type Description Default
T ndarray

Temperature grid :math:T in kelvin.

required
Im float

Peak intensity :math:I_m.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature :math:T_m in kelvin.

required
b float

Kinetic order parameter.

required

Returns:

Type Description
ndarray

Intensity values :math:I(T).

References

.. [1] Kitis, G., Gomes-Ros, J. M., and Tuyn, J. W. N. (1998). Thermoluminescence glow-curve deconvolution functions for first, second and general orders of kinetics.

Source code in tldecpy/models/go.py
def go_kg(T: np.ndarray, Im: float, E: float, Tm: float, b: float) -> np.ndarray:
    r"""
    Evaluate the Kitis-general expression for general-order TL kinetics.

    Parameters
    ----------
    T : numpy.ndarray
        Temperature grid :math:`T` in kelvin.
    Im : float
        Peak intensity :math:`I_m`.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature :math:`T_m` in kelvin.
    b : float
        Kinetic order parameter.

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.

    References
    ----------
    .. [1] Kitis, G., Gomes-Ros, J. M., and Tuyn, J. W. N. (1998).
           Thermoluminescence glow-curve deconvolution functions for first,
           second and general orders of kinetics.
    """
    if abs(b - 1.0) < 1e-3:
        b = 1.001

    temperature = np.maximum(np.asarray(T, dtype=float), 1e-5)
    k = KB_EV

    delta = (2.0 * k * temperature) / E
    delta_m = (2.0 * k * Tm) / E

    arg = (E / (k * temperature)) * ((temperature - Tm) / Tm)
    exp_arg = np.exp(np.clip(arg, -50.0, 50.0))

    prefactor = Im * (b ** (b / (b - 1.0))) * exp_arg
    bracket = (
        (b - 1.0) * (1.0 - delta) * (temperature**2 / Tm**2) * exp_arg
        + 1.0
        + (b - 1.0) * delta_m
    )
    bracket = np.maximum(bracket, 1e-9)
    return prefactor * (bracket ** (-b / (b - 1.0)))

tldecpy.models.go.go_rq

go_rq(T, Im, E, Tm, b)

Evaluate the rational-quadratic approximation for general-order TL kinetics.

Parameters:

Name Type Description Default
T ndarray

Temperature grid :math:T in kelvin.

required
Im float

Peak intensity :math:I_m.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature :math:T_m in kelvin.

required
b float

Kinetic order parameter.

required

Returns:

Type Description
ndarray

Intensity values :math:I(T).

References

.. [1] Gomez-Ros, J. M., and Kitis, G. (2002). Computerised glow-curve deconvolution using general and mixed-order kinetics.

Source code in tldecpy/models/go.py
def go_rq(T: np.ndarray, Im: float, E: float, Tm: float, b: float) -> np.ndarray:
    r"""
    Evaluate the rational-quadratic approximation for general-order TL kinetics.

    Parameters
    ----------
    T : numpy.ndarray
        Temperature grid :math:`T` in kelvin.
    Im : float
        Peak intensity :math:`I_m`.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature :math:`T_m` in kelvin.
    b : float
        Kinetic order parameter.

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.

    References
    ----------
    .. [1] Gomez-Ros, J. M., and Kitis, G. (2002).
           Computerised glow-curve deconvolution using general and mixed-order kinetics.
    """
    if abs(b - 1.0) < 1e-3:
        b = 1.001

    temperature = np.asarray(T, dtype=float)
    k = KB_EV

    def F_func(x: np.ndarray) -> np.ndarray:
        num = np.polyval([1.0, 2.334733, 0.250621], x)
        den = np.polyval([1.0, 3.330657, 1.681534], x)
        return 1.0 - (num / den)

    x_t = E / (k * temperature)
    x_tm = E / (k * Tm)
    x_tm_arr = np.asarray([x_tm], dtype=float)

    exp_diff = np.exp(np.clip(x_tm - x_t, -50.0, 50.0))
    term_inside = (temperature / Tm) * exp_diff * F_func(x_t) - F_func(x_tm_arr)[0]
    bracket = 1.0 + ((b - 1.0) / b) * x_tm * term_inside
    bracket = np.maximum(bracket, 1e-9)

    return Im * exp_diff * (bracket ** (-b / (b - 1.0)))

tldecpy.models.go.go_ge

go_ge(T, Im, E, Tm, b)

Evaluate the exponentialized approximation for general-order TL kinetics.

Parameters:

Name Type Description Default
T ndarray

Temperature grid :math:T in kelvin.

required
Im float

Peak intensity :math:I_m.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature :math:T_m in kelvin.

required
b float

Kinetic order parameter.

required

Returns:

Type Description
ndarray

Intensity values :math:I(T).

References

.. [1] Gomez-Ros, J. M., and Kitis, G. (2002). Computerised glow-curve deconvolution using general and mixed-order kinetics.

Source code in tldecpy/models/go.py
def go_ge(T: np.ndarray, Im: float, E: float, Tm: float, b: float) -> np.ndarray:
    r"""
    Evaluate the exponentialized approximation for general-order TL kinetics.

    Parameters
    ----------
    T : numpy.ndarray
        Temperature grid :math:`T` in kelvin.
    Im : float
        Peak intensity :math:`I_m`.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature :math:`T_m` in kelvin.
    b : float
        Kinetic order parameter.

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.

    References
    ----------
    .. [1] Gomez-Ros, J. M., and Kitis, G. (2002).
           Computerised glow-curve deconvolution using general and mixed-order kinetics.
    """
    if abs(b - 1.0) < 1e-3:
        b = 1.001

    temperature = np.asarray(T, dtype=float)
    k = KB_EV

    arg = (E / (k * Tm**2)) * (temperature - Tm)
    exp_arg = np.exp(np.clip(arg, -50.0, 50.0))

    numerator = Im * exp_arg
    denominator_term = (1.0 / b) + ((b - 1.0) / b) * exp_arg
    denominator_term = np.maximum(denominator_term, 1e-9)
    return numerator * (denominator_term ** (-b / (b - 1.0)))

Mixed-order (MO)

tldecpy.models.mixed.mo_kitis

mo_kitis(T, Im, E, Tm, alpha)

Mixed-order expression after Kitis and Gomez-Ros (2000).

Source code in tldecpy/models/mixed.py
def mo_kitis(T: np.ndarray, Im: float, E: float, Tm: float, alpha: float) -> np.ndarray:
    """
    Mixed-order expression after Kitis and Gomez-Ros (2000).
    """
    alpha = _check_alpha(alpha)
    k = KB_EV
    T = np.maximum(T, 1e-5)

    Am = 2.6 - 1.9203 * alpha + 0.324 * (alpha**3.338)
    Rm = (Am + alpha) / (Am - alpha)

    delta = (2 * k * T) / E
    delta_m = (2 * k * Tm) / E

    arg_exp = (E / (k * T)) * ((T - Tm) / Tm)
    exp_theta = np.exp(np.clip(arg_exp, -50, 50))

    term1 = (np.exp((1.0 - delta_m) / Rm) - alpha) ** 2
    term3_arg = (T**2 / (Tm**2 * Rm)) * exp_theta * (1.0 - delta)
    term3 = np.exp(np.clip(term3_arg, -50, 50))

    numerator = Im * term1 * exp_theta * term3
    denom_pre = np.exp((1.0 - delta_m) / Rm)
    denom_bracket = term3 - alpha

    denominator = denom_pre * (denom_bracket**2)
    return numerator / np.maximum(denominator, 1e-9)

tldecpy.models.mixed.mo_quad

mo_quad(T, Im, E, Tm, alpha)

Quadratic mixed-order approximation after Gomez-Ros and Kitis (2002). Corrected with auto-normalization to ensure peak height = Im.

Source code in tldecpy/models/mixed.py
def mo_quad(T: np.ndarray, Im: float, E: float, Tm: float, alpha: float) -> np.ndarray:
    """
    Quadratic mixed-order approximation after Gomez-Ros and Kitis (2002).
    Corrected with auto-normalization to ensure peak height = Im.
    """
    alpha = _check_alpha(alpha)
    k = KB_EV
    T = np.maximum(T, 1e-5)

    Rm = (1.0 - alpha) * (1.0 + 0.2922 * alpha - 0.2783 * alpha**2)
    z_Tm = E / (k * Tm)
    F_Tm = exp_int_approx(z_Tm)

    # Calculate shape function S(T)
    def calculate_shape(temp_arr):
        z = E / (k * temp_arr)
        F_val = exp_int_approx(z)

        # exp( E/kTm - E/kT )
        exp_diff = np.exp(np.clip(z_Tm - z, -50, 50))

        bracket_arg = Rm * z_Tm * ((temp_arr / Tm) * exp_diff * F_val - F_Tm)
        exp_bracket = np.exp(np.clip(bracket_arg, -50, 50))

        C = (1.0 - Rm) / (1.0 + Rm)

        num = 4.0 * (Rm**2) * exp_diff * exp_bracket
        den = (1.0 + Rm) * ((exp_bracket - C) ** 2)
        return num / np.maximum(den, 1e-9)

    S_T = calculate_shape(T)

    # Normalize by value at Tm
    S_Tm = calculate_shape(np.array([Tm]))[0]

    return Im * (S_T / S_Tm)

tldecpy.models.mixed.mo_vej

mo_vej(T, Im, E, Tm, alpha)

Mixed-order expression by Vejnovic et al. (2008).

Notes

This implementation follows the four-parameter fitting form shown as Eq. (20) in Vejnovic et al., Radiation Measurements 43 (2008) 1325-1330, using:

  • :math:l = 2/(2-\alpha) (Eq. 21)
  • :math:\sigma = T - T_m
  • :math:\Delta = 2kT/E
Source code in tldecpy/models/mixed.py
def mo_vej(T: np.ndarray, Im: float, E: float, Tm: float, alpha: float) -> np.ndarray:
    r"""
    Mixed-order expression by Vejnovic et al. (2008).

    Notes
    -----
    This implementation follows the four-parameter fitting form shown as Eq. (20)
    in Vejnovic et al., Radiation Measurements 43 (2008) 1325-1330, using:

    - :math:`l = 2/(2-\alpha)` (Eq. 21)
    - :math:`\sigma = T - T_m`
    - :math:`\Delta = 2kT/E`
    """
    alpha = _check_alpha(alpha)
    k = KB_EV
    T = np.maximum(T, 1e-5)

    l_val = 2.0 / (2.0 - alpha)
    delta_T = (2 * k * T) / E

    arg_rise = (E / (k * Tm)) * ((T - Tm) / Tm)
    exp_rise = np.exp(np.clip(arg_rise, -50, 50))

    factor_l = (2.0 / l_val) - 1.0
    bracket_arg = (T**2 / Tm**2) * factor_l * exp_rise * (1.0 - delta_T)
    exp_bracket = np.exp(np.clip(bracket_arg, -50, 50))

    prefactor = alpha * Im * ((2.0 - l_val) ** 2) / (l_val - 1.0)

    numerator = prefactor * exp_bracket * exp_rise
    denominator = (exp_bracket - alpha) ** 2

    return numerator / np.maximum(denominator, 1e-9)

OTOR

tldecpy.models.otor_lw.otor_lw

otor_lw(T, Im, E, Tm, R, beta=1.0)

OTOR model using the Lambert W semi-analytical solution.

Implements the One Trap One Recombination center (OTOR) kinetics using the analytical Lambert W expressions derived by Kitis & Vlachos (2013) and refined by Sadek et al. (2015). The peak maximum is forced to occur at :math:T_m via empirical correction factors applied to the frequency factor :math:s.

Two branches are handled automatically:

  • :math:R < 1 branch — principal Lambert W function :math:W_0.
  • :math:R > 1 branch — branch :math:W_{-1}.

Parameters:

Name Type Description Default
T ndarray

Temperature array in kelvin (1-D, float64). Values below 1e-12 K are clamped to prevent division by zero.

required
Im float

Peak maximum intensity :math:I_m in detector counts. Must be > 0.

required
E float

Trap activation energy :math:E in eV. Range: (0.1, 5.0) eV.

required
Tm float

Temperature at peak maximum :math:T_m in kelvin. Range: > 0 K.

required
R float

Retrapping ratio :math:R = A_n / A_m, where :math:A_n is the retrapping coefficient and :math:A_m the recombination coefficient (both in cm³/s).

  • :math:R \to 0 — negligible retrapping (approaches FO kinetics).
  • :math:R = 0.5 — equal retrapping and recombination.
  • :math:R \to 1 — dominant retrapping (approaches SO kinetics).

.. warning:: A numerical singularity exists near :math:R \approx 0.96 on the :math:R < 1 branch. The implementation clamps :math:R \leq 0.95 internally. Always set the upper bound to at most 0.90 in :class:~tldecpy.schemas.PeakSpec bounds.

required
beta float

Heating rate :math:\beta in K/s. Kept for API compatibility; the current semi-analytical form absorbs :math:\beta into the empirical correction of :math:s and does not use it directly.

1.0

Returns:

Type Description
ndarray

TL intensity array :math:I(T) in the same units as Im. All values are non-negative (clipped to 0).

Notes

The intensity is computed as:

.. math::

I(T) = I_m \,
\exp\!\left[-\frac{E}{kT}\cdot\frac{T_m-T}{T_m}\right]
\frac{W(T_m) + W(T_m)^2}{W(T) + W(T)^2}

where :math:W(T) is the Lambert W function evaluated at a temperature-dependent argument derived from the OTOR rate equations.

Examples:

>>> import numpy as np
>>> import tldecpy as tl
>>> T = np.linspace(350.0, 600.0, 300)
>>> I = tl.get_model("otor_lw")(T, Im=1000.0, E=1.3, Tm=460.0, R=0.15)
>>> print(f"Peak at T={T[I.argmax()]:.1f} K, I_max={I.max():.1f}")
References

.. [1] Kitis, G., and Vlachos, N. D. (2013). General semi-analytical expressions for TL, OSL and other luminescence stimulation modes derived from the OTOR model using the Lambert-W function. Radiat. Meas. 48, 47. .. [2] Sadek, A. M., et al. (2015). Characterization of the one trap one recombination (OTOR) level model using the Lambert-W function. Appl. Radiat. Isot. 95, 214.

Source code in tldecpy/models/otor_lw.py
def otor_lw(
    T: np.ndarray, Im: float, E: float, Tm: float, R: float, beta: float = 1.0
) -> np.ndarray:
    r"""
    OTOR model using the Lambert W semi-analytical solution.

    Implements the One Trap One Recombination center (OTOR) kinetics using
    the analytical Lambert W expressions derived by Kitis & Vlachos (2013)
    and refined by Sadek et al. (2015).  The peak maximum is forced to occur
    at :math:`T_m` via empirical correction factors applied to the frequency
    factor :math:`s`.

    Two branches are handled automatically:

    - :math:`R < 1` branch — principal Lambert W function :math:`W_0`.
    - :math:`R > 1` branch — branch :math:`W_{-1}`.

    Parameters
    ----------
    T : numpy.ndarray
        Temperature array in kelvin (1-D, float64).  Values below 1e-12 K
        are clamped to prevent division by zero.
    Im : float
        Peak maximum intensity :math:`I_m` in detector counts.
        Must be > 0.
    E : float
        Trap activation energy :math:`E` in eV.  Range: (0.1, 5.0) eV.
    Tm : float
        Temperature at peak maximum :math:`T_m` in kelvin.  Range: > 0 K.
    R : float
        Retrapping ratio :math:`R = A_n / A_m`, where :math:`A_n` is the
        retrapping coefficient and :math:`A_m` the recombination coefficient
        (both in cm³/s).

        - :math:`R \to 0` — negligible retrapping (approaches FO kinetics).
        - :math:`R = 0.5` — equal retrapping and recombination.
        - :math:`R \to 1` — dominant retrapping (approaches SO kinetics).

        .. warning::
           A numerical singularity exists near :math:`R \approx 0.96` on the
           :math:`R < 1` branch.  The implementation clamps :math:`R \leq 0.95`
           internally.  **Always set the upper bound to at most 0.90** in
           :class:`~tldecpy.schemas.PeakSpec` bounds.
    beta : float, default=1.0
        Heating rate :math:`\beta` in K/s.  Kept for API compatibility;
        the current semi-analytical form absorbs :math:`\beta` into the
        empirical correction of :math:`s` and does not use it directly.

    Returns
    -------
    numpy.ndarray
        TL intensity array :math:`I(T)` in the same units as ``Im``.
        All values are non-negative (clipped to 0).

    Notes
    -----
    The intensity is computed as:

    .. math::

        I(T) = I_m \,
        \exp\!\left[-\frac{E}{kT}\cdot\frac{T_m-T}{T_m}\right]
        \frac{W(T_m) + W(T_m)^2}{W(T) + W(T)^2}

    where :math:`W(T)` is the Lambert W function evaluated at a
    temperature-dependent argument derived from the OTOR rate equations.

    Examples
    --------
    >>> import numpy as np
    >>> import tldecpy as tl
    >>> T = np.linspace(350.0, 600.0, 300)
    >>> I = tl.get_model("otor_lw")(T, Im=1000.0, E=1.3, Tm=460.0, R=0.15)
    >>> print(f"Peak at T={T[I.argmax()]:.1f} K, I_max={I.max():.1f}")

    References
    ----------
    .. [1] Kitis, G., and Vlachos, N. D. (2013). *General semi-analytical
           expressions for TL, OSL and other luminescence stimulation modes
           derived from the OTOR model using the Lambert-W function.*
           Radiat. Meas. 48, 47.
    .. [2] Sadek, A. M., et al. (2015). *Characterization of the one trap
           one recombination (OTOR) level model using the Lambert-W function.*
           Appl. Radiat. Isot. 95, 214.
    """
    k = KB_EV

    T = np.asarray(T, dtype=np.float64)
    T = np.maximum(T, 1e-12)

    # ---- Integral function F(T, E) ----
    # F(T,E) = T*exp(-E/kT) + (E/k)*Ei(-E/kT)
    F_T = expi_integral(T, E, k)
    F_Tm = expi_integral(np.array([Tm], dtype=np.float64), E, k)[0]

    # ---- Branch handling ----
    # We decide branch based on the incoming R, but we still keep it away from exactly 1.
    R = clamp_R_near_one(float(R))

    if R < 1.0:
        # === R < 1 branch (Eq 10.1 style) ===
        # Keep away from the empirical singularity near ~0.96
        R = float(np.clip(R, 1e-12, 0.95))

        # Empirical correction factor C_R = 1 - 1.05 * R^1.26
        C_R = 1.0 - 1.05 * (R**1.26)
        C_R = max(C_R, 1e-12)

        # term_R = R/(1-R) - ln((1-R)/R)
        term_R = (R / (1.0 - R)) - np.log((1.0 - R) / R)

        # Q = E * exp(E/kTm) / (k * Tm^2 * C_R)  (computed in log-space)
        # log_Q = ln(E) + E/(kTm) - ln(k*Tm^2*C_R)
        # Use exp(clipped) to avoid overflow
        log_Q = np.log(max(E, 1e-300)) + (E / (k * Tm)) - np.log(max(k * (Tm**2) * C_R, 1e-300))
        Q = np.exp(np.clip(log_Q, -700.0, 700.0))

        Z_Tm = term_R + Q * F_Tm
        Z_T = term_R + Q * F_T

        # W(exp(Z)) on principal branch, computed stably
        W_Tm = _wexpz_principal_stable(np.array([Z_Tm], dtype=np.float64))[0]
        W_T = _wexpz_principal_stable(Z_T)

    else:
        # === R > 1 branch (Eq 10.2 style) ===
        # Keep away from exactly 1
        if R <= 1.0:
            R = 1.0 + 1e-4

        # Correction factor: 2.963 - 3.24 * R^(-0.74)
        corr = 2.963 - 3.24 * (R**-0.74)
        # Keep the original sign behavior, but avoid a near-zero denominator
        if abs(corr) < 1e-12:
            corr = 1e-12 * (1.0 if corr >= 0 else -1.0)

        # Keep the established Q form here (to avoid changing behavior near the R>1 correction sign),
        # but still avoid exp overflow via safe_exp.
        exp_term = float(safe_exp(np.asarray([E / (k * Tm)], dtype=float))[0])
        Q = (E * exp_term) / (k * (Tm**2) * corr)

        # Z0 = |R/(1-R)| + log(|(1-R)/R|)
        abs_frac = R / (R - 1.0)
        Z0 = abs_frac + np.log((R - 1.0) / R)

        Z_T = Z0 + Q * F_T
        Z_Tm = Z0 + Q * F_Tm

        # Argument for branch -1: -exp(-Z)
        arg_T = -safe_exp(-Z_T)
        arg_Tm = float(-safe_exp(np.asarray([-Z_Tm], dtype=float))[0])

        # Real-valued W_{-1}(x) exists for x in [-1/e, 0).
        # Clip to stay in-domain.
        eps = 1e-15
        lower = -1.0 / np.e + eps
        upper = -1e-300
        arg_T = np.clip(arg_T, lower, upper)
        arg_Tm = float(np.clip(arg_Tm, lower, upper))

        W_T = lambertw_safe(arg_T, k=-1)
        W_Tm = lambertw_safe(arg_Tm, k=-1)

    # ---- Intensity ----
    # exp_arg = (-E/(kT)) * ((Tm - T)/Tm) = -E/(kT) + E/(kTm)
    exp_arg = (-E / (k * T)) * ((Tm - T) / Tm)
    term_exp = np.exp(np.clip(exp_arg, -700.0, 700.0))

    num = W_Tm + W_Tm**2
    den = W_T + W_T**2
    den = np.maximum(den, 1e-12)

    I_calc = Im * term_exp * (num / den)
    return np.maximum(I_calc, 0.0)

tldecpy.models.otor_wo.otor_wo

otor_wo(T, Im, E, Tm, R, beta=1.0)

Evaluate OTOR intensity using Wright Omega representation.

Parameters:

Name Type Description Default
T ndarray

Temperature grid in kelvin.

required
Im float

Peak intensity.

required
E float

Activation energy in eV.

required
Tm float

Peak temperature in kelvin.

required
R float

Retrapping ratio :math:R=A_n/A_m.

required
beta float

Heating rate (kept for API compatibility).

1.0

Returns:

Type Description
ndarray

Intensity values :math:I(T).

Source code in tldecpy/models/otor_wo.py
def otor_wo(
    T: np.ndarray, Im: float, E: float, Tm: float, R: float, beta: float = 1.0
) -> np.ndarray:
    """
    Evaluate OTOR intensity using Wright Omega representation.

    Parameters
    ----------
    T : numpy.ndarray
        Temperature grid in kelvin.
    Im : float
        Peak intensity.
    E : float
        Activation energy in eV.
    Tm : float
        Peak temperature in kelvin.
    R : float
        Retrapping ratio :math:`R=A_n/A_m`.
    beta : float, default=1.0
        Heating rate (kept for API compatibility).

    Returns
    -------
    numpy.ndarray
        Intensity values :math:`I(T)`.
    """
    k = KB_EV
    R = clamp_R_near_one(R)

    # 1. Integral F(T)
    F_T = expi_integral(T, E, k)
    F_Tm = expi_integral(np.array([Tm]), E, k)[0]

    # 2. Correction factor (same empirical branch logic used by the LW variant)
    # Wright Omega form uses the same denominator correction
    # (1 - 1.05 * R^1.26) typically for R < 1.
    # Note: Wright Omega handles the transition continuously, but the empirical
    # correction factor for peak position derived by Kitis depends on R.
    # We use the same logic as otor_lw for the 'Q' factor.

    if R < 1.0:
        corr = 1.0 - 1.05 * (R**1.26)
    else:
        corr = 2.963 - 3.24 * (R**-0.74)

    if abs(corr) < 1e-3:
        corr = 1e-3 * np.sign(corr)

    exp_term = float(safe_exp(np.asarray([E / (k * Tm)], dtype=float))[0])
    Q = (E * exp_term) / (k * Tm**2 * corr)

    # 3. Argument Z for Wright Omega.
    # Wright Omega solves w + ln(w) = z and corresponds to W(exp(z)).

    if R < 1:
        # Term: ln( (1-R)/R )
        log_term = np.log((1.0 - R) / R)
        Z0 = (R / (1.0 - R)) - log_term
    else:
        # For R > 1, use the same real-valued branch convention as the LW variant.
        abs_frac = R / (R - 1.0)  # Negative number
        log_term = np.log((R - 1.0) / R)
        Z0 = abs_frac + log_term

    Z_T = Z0 + Q * F_T
    Z_Tm = Z0 + Q * F_Tm

    # 4. Evaluate Wright Omega.
    w_T = wrightomega_safe(Z_T)
    w_Tm = wrightomega_safe(Z_Tm)

    # 5. Intensity
    term_exp = safe_exp((-E / (k * T)) + (E / (k * Tm)))

    num = w_Tm + w_Tm**2
    den = w_T + w_T**2

    # Guard
    den = np.where(den == 0, 1e-9, den)

    I_calc = Im * term_exp * (num / den)

    return np.maximum(I_calc, 0.0)

Continuous trap distributions

tldecpy.models.continuous.continuous_gaussian

continuous_gaussian(T, Tn=None, In=None, E0=None, sigma=0.05, *, Im=None, E=None, Tm=None, n_energy=801, k=KB_EV)

Evaluate TL intensity for a Gaussian continuous trap-energy distribution.

Parameters:

Name Type Description Default
T ArrayLike

Temperature grid :math:T in kelvin.

required
Tn float | None

Characteristic temperature :math:T_N (K).

None
In float | None

Characteristic intensity :math:I_N.

None
E0 float | None

Characteristic activation energy :math:E_0 (eV).

None
sigma float

Standard deviation of the Gaussian energy distribution (eV).

0.05
Im float | None

Legacy aliases mapped internally to In, E0 and Tn.

None
E float | None

Legacy aliases mapped internally to In, E0 and Tn.

None
Tm float | None

Legacy aliases mapped internally to In, E0 and Tn.

None
n_energy int

Number of quadrature nodes for the energy integral.

801
k float

Boltzmann constant in eV/K.

KB_EV

Returns:

Type Description
ndarray

Simulated TL intensity evaluated at each value of T.

Source code in tldecpy/models/continuous.py
def continuous_gaussian(
    T: ArrayLike,
    Tn: float | None = None,
    In: float | None = None,
    E0: float | None = None,
    sigma: float = 0.05,
    *,
    Im: float | None = None,
    E: float | None = None,
    Tm: float | None = None,
    n_energy: int = 801,
    k: float = KB_EV,
) -> FloatArray:
    r"""
    Evaluate TL intensity for a Gaussian continuous trap-energy distribution.

    Parameters
    ----------
    T : ArrayLike
        Temperature grid :math:`T` in kelvin.
    Tn : float | None, optional
        Characteristic temperature :math:`T_N` (K).
    In : float | None, optional
        Characteristic intensity :math:`I_N`.
    E0 : float | None, optional
        Characteristic activation energy :math:`E_0` (eV).
    sigma : float, default=0.05
        Standard deviation of the Gaussian energy distribution (eV).
    Im, E, Tm : float | None, optional
        Legacy aliases mapped internally to ``In``, ``E0`` and ``Tn``.
    n_energy : int, default=801
        Number of quadrature nodes for the energy integral.
    k : float, default=KB_EV
        Boltzmann constant in eV/K.

    Returns
    -------
    numpy.ndarray
        Simulated TL intensity evaluated at each value of ``T``.
    """
    resolved_Tn, resolved_In, resolved_E0 = _resolve_continuous_aliases(
        Tn=Tn,
        In=In,
        E0=E0,
        Im=Im,
        E=E,
        Tm=Tm,
    )
    return continuous_glow_peak(
        T=T,
        Tn=resolved_Tn,
        In=resolved_In,
        E0=resolved_E0,
        sigma=sigma,
        model="continuous_gaussian",
        n_energy=n_energy,
        k=k,
    )

tldecpy.models.continuous.continuous_exponential

continuous_exponential(T, Tn=None, In=None, E0=None, sigma=0.05, *, Im=None, E=None, Tm=None, n_energy=801, k=KB_EV)

Evaluate TL intensity for an Exponential continuous trap-energy distribution.

Parameters:

Name Type Description Default
T ArrayLike

Temperature grid :math:T in kelvin.

required
Tn float | None

Characteristic temperature :math:T_N (K).

None
In float | None

Characteristic intensity :math:I_N.

None
E0 float | None

Characteristic activation energy :math:E_0 (eV).

None
sigma float

Characteristic width of the exponential tail (eV).

0.05
Im float | None

Legacy aliases mapped internally to In, E0 and Tn.

None
E float | None

Legacy aliases mapped internally to In, E0 and Tn.

None
Tm float | None

Legacy aliases mapped internally to In, E0 and Tn.

None
n_energy int

Number of quadrature nodes for the energy integral.

801
k float

Boltzmann constant in eV/K.

KB_EV

Returns:

Type Description
ndarray

Simulated TL intensity evaluated at each value of T.

Source code in tldecpy/models/continuous.py
def continuous_exponential(
    T: ArrayLike,
    Tn: float | None = None,
    In: float | None = None,
    E0: float | None = None,
    sigma: float = 0.05,
    *,
    Im: float | None = None,
    E: float | None = None,
    Tm: float | None = None,
    n_energy: int = 801,
    k: float = KB_EV,
) -> FloatArray:
    r"""
    Evaluate TL intensity for an Exponential continuous trap-energy distribution.

    Parameters
    ----------
    T : ArrayLike
        Temperature grid :math:`T` in kelvin.
    Tn : float | None, optional
        Characteristic temperature :math:`T_N` (K).
    In : float | None, optional
        Characteristic intensity :math:`I_N`.
    E0 : float | None, optional
        Characteristic activation energy :math:`E_0` (eV).
    sigma : float, default=0.05
        Characteristic width of the exponential tail (eV).
    Im, E, Tm : float | None, optional
        Legacy aliases mapped internally to ``In``, ``E0`` and ``Tn``.
    n_energy : int, default=801
        Number of quadrature nodes for the energy integral.
    k : float, default=KB_EV
        Boltzmann constant in eV/K.

    Returns
    -------
    numpy.ndarray
        Simulated TL intensity evaluated at each value of ``T``.
    """
    resolved_Tn, resolved_In, resolved_E0 = _resolve_continuous_aliases(
        Tn=Tn,
        In=In,
        E0=E0,
        Im=Im,
        E=E,
        Tm=Tm,
    )
    return continuous_glow_peak(
        T=T,
        Tn=resolved_Tn,
        In=resolved_In,
        E0=resolved_E0,
        sigma=sigma,
        model="continuous_exponential",
        n_energy=n_energy,
        k=k,
    )

Registry utilities

tldecpy.models.registry.get_model

get_model(key)

Resolve and return a callable model implementation.

Parameters:

Name Type Description Default
key str

Canonical key, alias, or family key.

required

Returns:

Type Description
Callable

Numerical model function that maps TL parameters to intensity.

Source code in tldecpy/models/registry.py
def get_model(key: str) -> ModelFunc:
    """
    Resolve and return a callable model implementation.

    Parameters
    ----------
    key : str
        Canonical key, alias, or family key.

    Returns
    -------
    collections.abc.Callable
        Numerical model function that maps TL parameters to intensity.
    """
    return get_info(key).func

tldecpy.models.registry.list_models

list_models(order=None, include_aliases=False)

List available model keys for TL kinetic families.

Parameters:

Name Type Description Default
order str | None

Optional family filter (fo, so, go, mo, otor, cont). Aliases (mix, continuous, dist) are also accepted.

None
include_aliases bool

When True, include legacy alias keys in the output.

False

Returns:

Type Description
list[str]

Canonical model keys (and optional aliases).

Source code in tldecpy/models/registry.py
def list_models(order: str | None = None, include_aliases: bool = False) -> List[str]:
    """
    List available model keys for TL kinetic families.

    Parameters
    ----------
    order : str | None, optional
        Optional family filter (``fo``, ``so``, ``go``, ``mo``, ``otor``, ``cont``).
        Aliases (``mix``, ``continuous``, ``dist``) are also accepted.
    include_aliases : bool, default=False
        When ``True``, include legacy alias keys in the output.

    Returns
    -------
    list[str]
        Canonical model keys (and optional aliases).
    """
    if order is None:
        keys = list(CANONICAL_MODELS.keys())
        if include_aliases:
            keys.extend(sorted(_ALIAS_TO_CANONICAL.keys()))
        return keys

    family_key = order
    if family_key in {"mix"}:
        family_key = "mo"
    if family_key in {"continuous", "dist"}:
        family_key = "cont"

    if family_key not in _FAMILY_MEMBERS:
        raise ValueError(f"Unknown order: {order}. Try: {list(_FAMILY_MEMBERS.keys())}")

    keys = list(_FAMILY_MEMBERS[family_key])
    if include_aliases:
        alias_keys = [
            alias
            for alias, canonical in _ALIAS_TO_CANONICAL.items()
            if CANONICAL_MODELS[canonical].family == family_key
        ]
        keys.extend(sorted(alias_keys))
    return keys

tldecpy.models.registry.get_model_info

get_model_info(key)

Return metadata for a canonical model key, alias, or family key.

Parameters:

Name Type Description Default
key str

Canonical key, alias, or family key (for example fo).

required

Returns:

Type Description
ModelInfo

Rich model metadata, including callable and display fields.

Source code in tldecpy/models/registry.py
def get_model_info(key: str) -> ModelInfo:
    """
    Return metadata for a canonical model key, alias, or family key.

    Parameters
    ----------
    key : str
        Canonical key, alias, or family key (for example ``fo``).

    Returns
    -------
    ModelInfo
        Rich model metadata, including callable and display fields.
    """
    return get_info(key)

tldecpy.models.registry.ModelInfo dataclass

ModelInfo(key, family, label, description, func, aliases=())

Metadata and callable for one model variant.