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

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

from typing import Any, Type
from warnings import warn

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

from .coupling import Coupling
from ..lib.utils import ArrayLike
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: float | Type[Coupling],
                     parameter: 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: tuple | list = None, keywords: dict = None):
        self.distribution = dist

        if parameter is not None:
            self.dist_parameter = parameter
        else:
            self.dist_parameter = ()

        if keywords is not None:
            self.dist_kw = keywords
        else:
            self.dist_kw = {}

    def t1(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
           mode: str = 'bpp', **kwargs) -> 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)
        else:
            return self.t1_csa(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) -> 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) -> 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) -> 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) -> 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(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
           mode: str = 'bpp', **kwargs) -> 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.t2_bpp(omega, tau, *specdens_args, **kwargs)
        elif mode == 'dipolar':
            return self.t2_dipolar(omega, tau, *specdens_args, **kwargs)
        else:
            return self.t2_csa(omega, tau, *specdens_args, **kwargs)

    def t2_bpp(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
               inverse: bool = True, prefactor: float = None) -> 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,
    ) -> 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,
    ) -> 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: Type[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 set_data(self, temp: ArrayLike, t1: ArrayLike):
        """
        Set data for evaluation.

        Args:
            temp (array-like): Temperature values
            t1 (array-like): T1 values

        """
        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: tuple | list = None,
            prefactor: tuple | list | float = None,
            coupling_kwargs: dict = None,
    ) -> None:
        """
        Determine a single parameter from a T1 minimum.
        It replaces the previously set value.

        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): Prefactor for
            coupling_kwargs (dict, optional): Keyword arguments for coupling, replacing old values

        Returns:
            A tuple of the value of varied parameter, or nan if nothing was varied
            and the minimum height calculated for given parameters.

        """

        stretching = minimon = 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
                prefactor = stretching

        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
            minimon = np.min(self.t1(omega, np.logspace(*tau_lims, num=1001), *self.dist_parameter,
                                     prefactor=self.prefactor))

        return stretching, minimon

    def calculate_t1_min(
            self,
            interpolate: int = None,
            trange: tuple[float, float] = None,
            use_log: bool = False,
    ) -> tuple[tuple[float, float], tuple[np.ndarray, np.ndarray] | None]:
        """
        Determine a minimum position for given T1 data

        Args:
            interpolate (int, optional):
              *  0 or None: No interpolation, minimum is data minimum
              *  1: Interpolation with a parabola
              *  2: Interpolation with a cubic spline
              *  3: Interpolation with Akima spline (less wiggly than cubic)
            trange (tuple): Range (left border, range border) of interpolation in K.
                Interpolation without a given range uses two points left and right of minimum value.
            use_log (bool):  Default is `True`.

        Returns:
            The minimum position (`T_min`, `T1_min`)

        """
        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: list | tuple = None,
            prefactor: float = None,
            coupling_param: list = None,
            coupling_kwargs: dict = None,
    ) -> tuple[np.ndarray, dict]:
        """
        Calculate correlation times from set T1 data.
        Optional arguments overwrite previousliy set parameter.

        Args:
            mode (str, {`raw`, `mean`, `logmean`, `max`}): Type of correlation time. Default is `raw`.
            interpolate (bool): If ``True`` and T1 minimum was determined by nterpolation,
                T1 on interpolated line instead of measured value is used. Default is `False`.
            omega (float, optional): Larmor frequency (in 1/s)
            dist_parameter (list, optional): List of parameter of spectral density
            prefactor (float, optional): Prefactor of T1 calculation, will
            coupling_param (list, optional): Parameter for coupling constant, ignored if `prefactor` is given.
            coupling_kwargs (dict, optional): Keyword arguments for coupling constant, ignored if `prefactor` is given.

        Returns:
            Array of temperature and correlation times and a dictionary of used parameter during calculation.

        """
        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
                coupling_param = self.coup_parameter
            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=1001)
        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(f'Value {t1_i} below set minimum, wonky correlation time')
                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=1001)
        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'Value {t1_i} below set minimum, wonky correlation time')
                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,
                'prefactor': self.prefactor}
        if self.coupling is not None:
            opts['coupling'] = (self.coupling.name, coupling_param, coupling_kwargs)

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