"""
Relaxation
==========

Classes to calculate spin-lattice and spin-spin relaxation, as well as to evaluate T1 data and calculate correlation times
"""

from pathlib import Path
from typing import Any, Optional, Tuple, Type, Union
from warnings import warn

import numpy as np
from scipy.interpolate import interp1d, Akima1DInterpolator
from scipy.optimize import minimize

from nmreval.lib.utils import ArrayLike
from .coupling import Coupling
from ..distributions.base import Distribution
from ..distributions.debye import Debye as Debye
from ..utils import gamma_full


class Relaxation:
    """
    Class to calculate relaxation times.
    """

    def __init__(self, distribution: Type[Distribution] = None):
        self.distribution = distribution
        self.dist_parameter = ()
        self._dist_kw = {}

        self.coupling = None
        self.coup_parameter = []
        self.coup_kw = {}

        self.prefactor = 1.

    def __repr__(self):
        if self.distribution is not None:
            return str(self.distribution.name)
        else:
            return super().__repr__()

    def set_coupling(self, coupling: Union[float, Type[Coupling]],
                     parameter: Union[tuple, list] = None, keywords: dict = None):

        if parameter is not None:
            self.coup_parameter = parameter

        if keywords is not None:
            self.coup_kw = keywords

        if isinstance(coupling, float):
            self.prefactor = coupling
        elif issubclass(coupling, Coupling):
            self.coupling = coupling
            if self.coup_kw is None:
                raise ValueError('Coupling is missing parameter')
            self.prefactor = self.coupling.relax(*self.coup_parameter, **self.coup_kw)
        else:
            raise ValueError(f'`coupling` is not number or of type `Coupling`, found {coupling!r}')

    def set_distribution(self, dist: Type[Distribution], parameter: Union[tuple, list] = None, keywords: dict = None):
        self.distribution = dist

        if parameter is not None:
            self.dist_parameter = parameter

        if keywords is not None:
            self._dist_kw = keywords

    def t1(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
           mode: str = 'bpp', **kwargs) -> Union[np.ndarray, float]:
        r"""
        Convenience function

        Args:
            omega (array-like): Frequency in 1/s (not Hz)
            tau (array-like): Correlation times in s
            *specdens_args: Additional parameter for spectral density.
                If given this will replace previously set arguments during calculation
            mode (str, {`bpp`, `dipolar`, `csa`}): Type of relaxation

                - `bpp` : Homonuclear dipolar/quadrupolar interaction :func:`(see here) <nmreval.nmr.Relaxation.t1_bpp>`
                - `dipolar` : Heteronuclear dipolar interaction :func:`(see here) <nmreval.nmr.Relaxation.t1_dipolar>`
                - `csa` : Chemical shift interaction :func:`(see here) <nmreval.nmr.Relaxation.t1_csa>`

                Default is `bpp`

            **kwargs:

        Returns:
            float or ndarray:
            A k by n array where k is length of ``omega`` and n is length of ``tau``.
            If both are of length 1 a number is returned instead of array.

        """
        if mode not in ['bpp', 'dipolar', 'csa']:
            raise ValueError(f'Unknown mode {mode} not `bpp`, `dipolar`, `csa`.')

        # defaults to BPP
        if mode == 'bpp':
            return self.t1_bpp(omega, tau, *specdens_args, **kwargs)
        elif mode == 'dipolar':
            return self.t1_dipolar(omega, tau, *specdens_args, **kwargs)

    def t1_dipolar(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any, inverse: bool = True,
                   prefactor: float = None, omega_coup: ArrayLike = None,
                   gamma_coup: str = None, gamma_obs: str = None) -> Union[np.ndarray, float]:
        r"""Calculate T1 under heteronuclear dipolar coupling.

        .. math::
            \frac{1}{T_1} = C [3J(\omega_I) + 6J(\omega_I + \omega_I) + J(\omega_I - \omega_S)]

        Args:
            omega (array-like): Frequency in 1/s (not Hz)
            tau (array-like): Correlation times in s
            *specdens_args: Additional parameter for spectral density.
                If given this will replace previously set arguments during calculation
            inverse (bool): Function returs relaxation time if True else relaxation rate. Default is `True`
            prefactor (float, optional): If given it is used as prefactor `C`
                instead of previously set coupling prefactor.
            omega_coup (array-like, optional): Frequency of coupled isotope must be of same shape as omega.
            gamma_obs (str, optional): Name of observed  isotope ('1H', '23Na', ...).
                Values is used if `omega_coup` is None.
            gamma_coup (str, optional): Name of coupled isotope ('1H', '23Na', ...).
                Values is used if `omega_coup` is None.

        Returns:
            float or ndarray:
            A k by n array where k is length of ``omega`` and n is length of ``tau``.
            If both are of length 1 a number is returned instead of array.

        """
        omega = np.asanyarray(omega)

        if omega_coup is not None:
            omega_coup = np.asanyarray(omega_coup)
            if omega.shape != omega_coup.shape:
                raise ValueError('omega and omega_coup have not same shape')
            omega_2 = omega_coup
        elif (gamma_coup is not None) and (gamma_obs is not None):
            omega_2 = gamma_full[gamma_coup] / gamma_full[gamma_obs] * omega
        else:
            raise AttributeError('Unknown second frequency. Set `omega_coup`, or `gamma_coup` and `gamma_obs`')

        if prefactor is None:
            prefactor = self.prefactor

        if len(specdens_args) == 0:
            specdens_args = self.dist_parameter

        rate = prefactor * (self.distribution.specdens(omega - omega_2, tau, *specdens_args) +
                            3 * self.distribution.specdens(omega, tau, *specdens_args) +
                            6 * self.distribution.specdens(omega + omega_2, tau, *specdens_args)) / 2
        if inverse:
            return 1 / rate
        else:
            return rate

    def t1_bpp(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
               inverse: bool = True, prefactor: float = None) -> Union[np.ndarray, float]:
        r"""Calculate T1 under homonuclear dipolar coupling or quadrupolar coupling.

        .. math::
            \frac{1}{T_1} = C [J(\omega) + 4J(2\omega)]

        Args:
            omega (array-like): Frequency in 1/s (not Hz)
            tau (array-like): Correlation times in s
            *specdens_args: Additional parameter for spectral density.
                If given this will replace previously set arguments during calculation
            inverse (bool): Function returs relaxation time if True else relaxation rate. Default is `True`
            prefactor (float, optional): If given it is used as prefactor `C`
                instead of previously set coupling prefactor.

        Returns:
            float or ndarray:
            A k by n array where k is length of ``omega`` and n is length of ``tau``.
            If both are of length 1 a number is returned instead of array.

        """
        if prefactor is None:
            prefactor = self.prefactor

        if len(specdens_args) == 0:
            specdens_args = self.dist_parameter

        rate = prefactor * (self.distribution.specdens(omega, tau, *specdens_args) +
                            4 * self.distribution.specdens(2 * omega, tau, *specdens_args))
        if inverse:
            return 1. / rate
        else:
            return rate

    def t1_csa(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
               inverse: bool = True, prefactor: float = None) -> Union[np.ndarray, float]:
        r"""Calculate T1 under chemical shift anisotropy.

        .. math::
            \frac{1}{T_1} = C \omega^2 [J(\omega)]

        This relation disregards antisymmetric CSA contributions

        Args:
            omega (array-like): Frequency in 1/s (not Hz)
            tau (array-like): Correlation times in s
            *specdens_args: Additional parameter for spectral density.
                If given this will replace previously set arguments during calculation
            inverse (bool): Function returs relaxation time if True else relaxation rate. Default is `True`
            prefactor (float, optional): If given it is used as prefactor `C`
                instead of previously set coupling prefactor.

        Returns:
            float or ndarray:
            A k by n array where k is length of ``omega`` and n is length of ``tau``.
            If both are of length 1 a number is returned instead of array.

        """
        if prefactor is None:
            prefactor = self.prefactor

        if len(specdens_args) == 0:
            specdens_args = self.dist_parameter

        rate = prefactor * (self.distribution.specdens(omega, tau, *specdens_args) +
                            4 * self.distribution.specdens(2 * omega, tau, *specdens_args))
        if inverse:
            return 1. / rate
        else:
            return rate

    def t1q(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
            inverse: bool = True, prefactor: float = None) -> Union[np.ndarray, float]:
        r"""Calculate T1q for homonuclear dipolar coupling or quadrupolar coupling (I=1).

        .. math::
            \frac{1}{T_{1q}} = 3 C J(\omega)

        Args:
            omega (array-like): Frequency in 1/s (not Hz)
            tau (array-like): Correlation times in s
            *specdens_args: Additional parameter for spectral density.
                If given this will replace previously set arguments during calculation
            inverse (bool): Function returs relaxation time if True else relaxation rate. Default is `True`
            prefactor (float, optional): If given it is used as prefactor `C`
                instead of previously set coupling prefactor.

        Returns:
            float or ndarray:
            A k by n array where k is length of ``omega`` and n is length of ``tau``.
            If both are of length 1 a number is returned instead of array.

        """
        if prefactor is None:
            prefactor = self.prefactor

        if len(specdens_args) == 0:
            specdens_args = self.dist_parameter

        rate = 3 * prefactor * self.distribution.specdens(omega, tau, *specdens_args)
        if inverse:
            return 1. / rate
        else:
            return rate

    def t1_rho(self, omega: ArrayLike, tau: ArrayLike, omega_rf: float, *specdens_args: Any,
               inverse: bool = True, prefactor: float = None):
        r"""Calculate T1rho for homonuclear dipolar coupling or quadrupolar coupling.

        .. math::
            \frac{1}{T_{1\rho}} = \frac{C}{2} [5J(\omega) + 2J(2\omega) + 3J(2\omega_{RF})]

        Args:
            omega (array-like): Frequency in 1/s (not Hz).
            tau (array-like): Correlation times in s.
            omega_rf (float): Frequency of RF lock IN 1/s.
            *specdens_args: Additional parameter for spectral density.
                If given this will replace previously set arguments during calculation
            inverse (bool): Function returs relaxation time if True else relaxation rate. Default is `True`
            prefactor (float, optional): If given it is used as prefactor `C`
                instead of previously set coupling prefactor.

        Returns:
            float or ndarray:
            A k by n array where k is length of ``omega`` and n is length of ``tau``.
            If both are of length 1 a number is returned instead of array.

        """
        if prefactor is None:
            prefactor = self.prefactor

        if len(specdens_args) == 0:
            specdens_args = self.dist_parameter

        rate = prefactor * (5 * self.distribution.specdens(omega, tau, *specdens_args) +
                            2 * self.distribution.specdens(2 * omega, tau, *specdens_args) +
                            3 * self.distribution.specdens(2 * omega_rf, tau, *specdens_args)) / 2
        if inverse:
            return 1. / rate
        else:
            return rate

    def t2_bpp(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
               inverse: bool = True, prefactor: float = None) -> Union[np.ndarray, float]:
        r"""Calculate T2 under homonuclear dipolar coupling or quadrupolar coupling.

        .. math::
            \frac{1}{T_2} = \frac{C}{2} [3J() + 5J(2\omega) + 3J(2\omega)]

        Args:
            omega (array-like): Frequency in 1/s (not Hz)
            tau (array-like): Correlation times in s
            *specdens_args (optional): Additional parameter for spectral density.
                If given this will replace previously set arguments during calculation
            inverse (bool): Function returs relaxation time if True else relaxation rate. Default is `True`
            prefactor (float, optional): If given it is used as prefactor `C`
                instead of previously set coupling prefactor.

        Returns:
            float or ndarray:
            A k by n array where k is length of ``omega`` and n is length of ``tau``.
            If both are of length 1 a number is returned instead of array.

        """
        if prefactor is None:
            prefactor = self.prefactor

        if len(specdens_args) == 0:
            specdens_args = self.dist_parameter

        rate = prefactor * (3 * self.distribution.specdens(0, tau, *specdens_args) +
                            5 * self.distribution.specdens(omega, tau, *specdens_args) +
                            2 * self.distribution.specdens(2 * omega, tau, *specdens_args)) / 2
        if inverse:
            return 1. / rate
        else:
            return rate

    def t2_dipolar(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
                   inverse: bool = True, prefactor: float = None, omega_coup: ArrayLike = None,
                   gamma_coup: str = None, gamma_obs: str = None) -> Union[np.ndarray, float]:
        r"""Calculate T2 under heteronuclear dipolar coupling.

        .. math::
            \frac{1}{T_2} = \frac{C}{2} [4J(0) + J(\omega_I - \omega_I) + 3J(\omega_S) + 6J(\omega_I) + 6J(\omega_I + \omega_I)]

        Args:
            omega (array-like): Frequency in 1/s (not Hz)
            tau (array-like): Correlation times in s
            *specdens_args (optional): Additional parameter for spectral density.
                If given this will replace previously set arguments during calculation
            inverse (bool): Function returs relaxation time if True else relaxation rate. Default is `True`
            prefactor (float, optional): If given it is used as prefactor `C`
                instead of previously set coupling prefactor.
            omega_coup (array-like, optional): Frequency of coupled isotope must be of same shape as omega.
            gamma_obs (str, optional): Name of observed  isotope ('1H', '23Na', ...).
                Values is used if `omega_coup` is None.
            gamma_coup (str, optional): Name of coupled isotope ('1H', '23Na', ...).
                Values is used if `omega_coup` is None.

        Returns:
            float or ndarray:
            A k by n array where k is length of ``omega`` and n is length of ``tau``.
            If both are of length 1 a number is returned instead of array.

        """
        omega = np.asanyarray(omega)

        if omega_coup is not None:
            omega_coup = np.asanyarray(omega_coup)
            if omega.shape != omega_coup.shape:
                raise ValueError('omega and omega_coup have not same shape')
            omega_2 = omega_coup
        elif (gamma_coup is not None) and (gamma_obs is not None):
            omega_2 = gamma_full[gamma_coup] / gamma_full[gamma_obs] * omega
        else:
            raise AttributeError('Unknown second frequency. Set `omega_coup`, or `gamma_coup` and `gamma_obs`')

        if prefactor is None:
            prefactor = self.prefactor

        if len(specdens_args) == 0:
            specdens_args = self.dist_parameter

        rate = prefactor * (4 * self.distribution.specdens(0, tau, *specdens_args) +
                            self.distribution.specdens(omega - omega_2, tau, *specdens_args) +
                            3 * self.distribution.specdens(omega_2, tau, *specdens_args) +
                            6 * self.distribution.specdens(omega, tau, *specdens_args) +
                            6 * self.distribution.specdens(omega + omega_2, tau, *specdens_args)) / 2.

        if inverse:
            return 1. / rate
        else:
            return rate

    def t2_csa(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
               inverse: bool = True, prefactor: float = None) -> Union[np.ndarray, float]:
        r"""Calculate T1 under chemical shift anisotropy.

        .. math::
            \frac{1}{T_2} = \frac{C}{2} \omega^2 \bigl[J(0) + \frac{4}{3}J(\omega)\bigr]

        Args:
            omega (array-like): Frequency in 1/s (not Hz)
            tau (array-like): Correlation times in s
            *specdens_args (optional): Additional parameter for spectral density.
                If given this will replace previously set arguments during calculation
            inverse (bool): Function returs relaxation time if True else relaxation rate. Default is `True`
            prefactor (float, optional): If given it is used as prefactor `C`
                instead of previously set coupling prefactor.

        Returns:
            float or ndarray:
            A k by n array where k is length of ``omega`` and n is length of ``tau``.
            If both are of length 1 a number is returned instead of array.

        """
        if prefactor is None:
            prefactor = self.prefactor

        if len(specdens_args) == 0:
            specdens_args = self.dist_parameter

        rate = prefactor * (self.distribution.specdens(0, tau, *specdens_args) +
                            4 / 3 * self.distribution.specdens(omega, tau, *specdens_args)) / 2.
        if inverse:
            return 1. / rate
        else:
            return rate


class RelaxationEvaluation(Relaxation):
    def __init__(self, distribution=None):
        super().__init__(distribution=distribution)
        self.t1min = (np.nan, np.nan)
        self._interpolate = None
        self._interpolate_range = (None, None)

        self.omega = np.nan
        self.x = None
        self.y = None

    def data(self, temp, t1):
        temp = np.asanyarray(temp)
        t1 = np.asanyarray(t1)
        sortidx = temp.argsort()
        self.x = temp[sortidx]
        self.y = t1[sortidx]
        self.calculate_t1_min()

    def get_increase(self, height: float = None, idx: int = 0, mode: str = None, omega: float = None,
                     dist_parameter: Union[tuple, list] = None, prefactor: Union[tuple, list, float] = None,
                     coupling_kwargs: dict = None):
        """
        Determine a single parameter from a T1 minimum

        Args:
            height (float, optional): Height of T1 minimum
            mode (str, {`distribution`, `prefactor`}, optional): If given,
            idx (int): Default is 0.
            omega (float, optional): Larmor frequency (in 1/s)
            dist_parameter (tuple, optional):
            prefactor (tuple, float, optional):
            coupling_kwargs (dict, optional):

        Returns:

        """

        stretching = mini = np.nan
        if height is None:
            height = self.t1min[1]

        if omega is None:
            if np.isnan(self.omega):
                raise ValueError('No larmor frequency set')

            omega = self.omega

        if prefactor is None:
            prefactor = self.prefactor

        if dist_parameter is None:
            dist_parameter = self.dist_parameter

        if coupling_kwargs is None:
            coupling_kwargs = self.coup_kw

        tau_lims = np.log10(1 / omega) - 3, np.log10(1 / omega) + 3

        if mode is None:
            # nothing is variable -> just calculate minimum for given parameter
            if isinstance(prefactor, (tuple, list)):
                if coupling_kwargs is None:
                    coupling_kwargs = self.coup_kw
                prefactor = self.coupling.relax(*prefactor, **coupling_kwargs)

            return stretching, np.min(self.t1(omega, np.logspace(*tau_lims, num=1001),
                                              *dist_parameter, prefactor=prefactor, inverse=True))

        if mode == 'distribution':
            # width parameter of spectral density is variable

            if isinstance(self.distribution, Debye):
                # Debye is easy
                return 1., np.min(self.t1(omega, np.logspace(*tau_lims, num=1001),
                                          *dist_parameter, prefactor=prefactor, inverse=True))

            if isinstance(prefactor, (list, tuple)):
                # calculate prefactor from coupling
                if self.coupling is None:
                    raise ValueError('Coupling must be set before evaluation of T1 min')

                if coupling_kwargs is None:
                    coupling_kwargs = self.coup_kw

                prefactor = self.coupling.relax(*prefactor, **coupling_kwargs)

            use_fmin = True
            if self.distribution.name in ['KWW', 'Cole-Davidson', 'Cole-Cole', 'Log-Gaussian']:
                # using precalculated values is faster than actual calculation, especially KWW and LG
                from importlib.resources import path
                with path('resources.nmr', 'T1_min.npz') as fp:
                    if fp.exists():
                        t1min_table = np.load(str(fp))[self.distribution.name.replace('-', '_')]

                        debye_min = omega / 1.42517571908650 / prefactor  # Zahl kommt von Wolfram alpha
                        ratio = height / debye_min
                        stretching = t1min_table[np.argmin(np.abs(t1min_table[:, 1] - ratio)), 0]

                        use_fmin = False
                        dist_parameter = [stretching]

            if use_fmin:
                # use for untabulated spectral densities or if something went wrong
                #
                def _f(_x, p, ii):
                    p[ii] = _x[0]
                    t1_calc = np.min(self.t1(omega, np.logspace(*tau_lims, num=1001), *p,
                                             prefactor=prefactor, inverse=True))

                    return np.abs(t1_calc - height)

                stretching = minimize(_f, np.array([dist_parameter[idx]]),
                                      args=(dist_parameter, idx),
                                      method='Nelder-Mead').x[0]
                dist_parameter[idx] = stretching

        elif mode == 'coupling':
            # variable is somewhere in the prefcotor

            t1_no_coup = np.min(self.t1(omega, np.logspace(*tau_lims, num=1001),
                                        *dist_parameter, prefactor=1))

            if isinstance(prefactor, (tuple, list)):
                prefactor = list(prefactor)

                def _f(_x, p, ii):
                    p[ii] = _x[0]
                    return np.abs(t1_no_coup / height - self.coupling.relax(*p, **coupling_kwargs))

                stretching = minimize(_f, np.array([prefactor[idx]]),
                                      args=(prefactor, idx)).x[0]

            else:
                stretching = t1_no_coup / height

        else:
            raise ValueError('Use `distribution` or `coupling` to set parameter')

        if stretching:
            self.prefactor = prefactor
            self.dist_parameter = dist_parameter
            if isinstance(prefactor, (tuple, list)):
                self.coup_parameter = prefactor
                self.prefactor = self.coupling.relax(*self.coup_parameter, **self.coup_kw)
            else:
                self.prefactor = prefactor
            mini = np.min(self.t1(omega, np.logspace(*tau_lims, num=1001), *self.dist_parameter,
                                  prefactor=self.prefactor))

        return stretching, mini

    def calculate_t1_min(self, interpolate: int = None, trange: Tuple[float, float] = None, use_log: bool = False) -> \
            Tuple[Tuple[float, float], Optional[Tuple[np.ndarray, np.ndarray]]]:
        min_index = np.argmin(self.y)
        t1_min = (self.x[min_index], self.y[min_index])
        parabola = None
        self._interpolate_range = (None, None)

        if interpolate is not None:
            if interpolate not in [0, 1, 2, 3]:
                raise ValueError(f'Unknown interpolation value {interpolate}')

            if interpolate != 0:
                if trange is None:
                    left_b = max(min_index - 2, 0)
                    right_b = min(len(self.x), min_index + 3)
                else:
                    left_b = np.argmin(np.abs(self.x - trange[0]))
                    right_b = np.argmin(np.abs(self.x - trange[1])) + 1

                temp_range = self.x[left_b:right_b]
                t1_range = self.y[left_b:right_b]
                if use_log:
                    t1_range = np.log10(t1_range)

                try:
                    x_inter = np.linspace(temp_range[0], temp_range[-1], num=51)

                    if interpolate == 1:
                        i_func = np.poly1d(np.polyfit(temp_range, t1_range, 2))
                    elif interpolate == 2:
                        i_func = interp1d(temp_range, t1_range, kind='cubic')
                    else:
                        i_func = Akima1DInterpolator(temp_range, t1_range)

                    if use_log:
                        y_inter = 10**i_func(x_inter)
                    else:
                        y_inter = i_func(x_inter)
                    t1_min = (x_inter[np.argmin(y_inter)], y_inter[np.argmin(y_inter)])
                    self._interpolate = i_func
                    parabola = (x_inter, y_inter)
                    self._interpolate_range = (temp_range[0], temp_range[-1])
                    t1_min = x_inter[np.argmin(y_inter)], y_inter[np.argmin(y_inter)]

                except IndexError:
                    pass

        self.t1min = t1_min

        return t1_min, parabola

    def correlation_from_t1(self, mode: str = 'raw', interpolate: bool = False, omega: float = None,
                            dist_parameter: Union[float, list, tuple] = None, prefactor: Union[float, tuple, list] = None,
                            coupling_param: list = None, coupling_kwargs: dict = None) -> Tuple[np.ndarray, dict]:

        if self.x is None:
            raise ValueError('Temperature is not set')

        if self.y is None:
            raise ValueError('T1 data is not set')

        if omega is None:
            if np.isnan(self.omega):
                raise ValueError('No frequency set')
            omega = self.omega

        if prefactor is None:
            if coupling_kwargs is None:
                coupling_kwargs = self.coup_kw

            if coupling_param is None:
                prefactor = self.prefactor
            else:
                prefactor = self.coupling.relax(*coupling_param, **coupling_kwargs)

        if dist_parameter is None:
            dist_parameter = self.dist_parameter

        fast_idx = self.x > self.t1min[0]

        slow_t1 = self.y[~fast_idx]
        slow_temp = self.x[~fast_idx]

        correlation_times = np.ones_like(self.x)
        offset = len(slow_t1)

        base_taus = np.logspace(-10, -7, num=1001)
        min_tau = base_taus[np.argmin(self.t1(omega, base_taus, *dist_parameter, prefactor=prefactor))]

        taus = np.geomspace(min_tau, 100. * min_tau, num=501)
        current_t1 = self.t1(omega, taus, *dist_parameter, prefactor=prefactor)

        for i in range(1, len(slow_t1) + 1):
            t1_i = slow_t1[-i]

            if interpolate and self._interpolate is not None:
                if slow_temp[-i] >= self._interpolate_range[0]:
                    t1_i = self._interpolate(slow_temp[-i])

            if np.min(current_t1) > t1_i:
                warn('Correlation time could not be calculated')
                correlation_times[offset - i] = taus[0]
                continue

            cross_idx = np.where(np.diff(np.sign(current_t1 - t1_i)))[0]
            while len(cross_idx) == 0:
                taus *= 10
                current_t1 = self.t1(omega, taus, *dist_parameter, prefactor=prefactor)
                cross_idx = np.where(np.diff(np.sign(current_t1 - t1_i)))[0]

            lamb = (t1_i - current_t1[cross_idx]) / (current_t1[cross_idx + 1] - current_t1[cross_idx])
            correlation_times[offset - i] = (taus[cross_idx + 1] * lamb + (1 - lamb) * taus[cross_idx])[0]

        fast_t1 = self.y[fast_idx]
        fast_temp = self.x[fast_idx]

        taus = np.geomspace(0.01 * min_tau, min_tau, num=501)
        current_t1 = self.t1(omega, taus, *dist_parameter, prefactor=prefactor)

        for i in range(len(fast_t1)):
            t1_i = fast_t1[i]

            if interpolate and self._interpolate is not None:
                if fast_temp[i] <= self._interpolate_range[1]:
                    t1_i = self._interpolate(fast_temp[i])

            if current_t1[-1] > t1_i:
                correlation_times[offset + i] = taus[-1]
                warn(f'Correlation time for {correlation_times[offset + i]} could not be calculated')
                continue

            cross_idx = np.where(np.diff(np.sign(current_t1 - t1_i)))[0]
            while len(cross_idx) == 0:
                taus /= 10
                current_t1 = self.t1(omega, taus, *dist_parameter, prefactor=prefactor)
                cross_idx = np.where(np.diff(np.sign(current_t1 - t1_i)))[0]

            lamb = (t1_i - current_t1[cross_idx]) / (current_t1[cross_idx + 1] - current_t1[cross_idx])
            correlation_times[offset + i] = (taus[cross_idx + 1] * lamb + (1 - lamb) * taus[cross_idx])[0]

        opts = {'distribution': (self.distribution.name, dist_parameter),
                'frequency': omega / 2 / np.pi}
        if self.coupling is not None:
            opts['coupling'] = (self.coupling.name, self.prefactor, coupling_param, coupling_kwargs)
        else:
            opts['coupling'] = (self.prefactor,)

        return np.c_[self.x, self.distribution.mean_value(correlation_times, *dist_parameter, mode=mode)], opts

    def _create_minimum_file(self, filename):
        x = np.geomspace(0.1, 1, num=10001)
        with Path(filename).open('w') as f:
            f.write('# broadening\tT1_min/min_Debye\n')
            for i, a in enumerate(np.geomspace(0.1, 10, num=10000)):
                alpha = a
                t1_i = self.t1_bpp(1, x, alpha)
                f.write(f'{alpha}\t{min(t1_i) / 0.701667864144}\n')
                f.flush()