import numpy as np

from scipy import integrate

from saf.euler1d.common.zndsolver import ZNDSolver

FORMAT = '22.16e'


class ScaleConverter:
    """Convert from Erpenbeck scales to scales of Sharpe and Lee&Stewart.

    Additionally, convert time quantities to Abouseif&Toong scales.

    Parameters
    ----------
    c : saf.euler1d.linear.config.Config
        Config object containing values in Erpenbeck scales.
    trace : bool (default: False)
        If True, the conversion factors will be printed to the standard output.

    """
    def __init__(self, c, trace=False):
        self.c = c
        self.trace = trace

        self._znd = ZNDSolver(self.c)

        L = np.ceil(self._znd.reaction_length)
        N = c.n12 * L
        x = np.linspace(-L, 0, num=N+1)
        self._znd.compute(x)

    def to_sharpe_scales(self):
        gamma, p_a, rho_a = self.c.gamma, self.c.p_a, self.c.rho_a
        d = self._znd.d

        gamma_sqrt = np.sqrt(gamma)
        c_a_sq = gamma * p_a / rho_a
        c_a = np.sqrt(c_a_sq)

        factors = {
            'density': 1.0,
            'volume': 1.0,
            'pressure': c_a_sq / (gamma * d**2),
            'velocity': c_a / (gamma_sqrt * d),
            'energy': c_a_sq / (gamma * d**2),
            'time': gamma_sqrt * d / c_a,
            'inv_time': c_a / (gamma_sqrt * d),
        }

        self._print_factors(factors, 'Sharpe')
        return factors

    def to_lee_and_stewart_scales(self):
        gamma, p_a, rho_a = self.c.gamma, self.c.p_a, self.c.rho_a

        gamma_sqrt = np.sqrt(gamma)
        c_a_sq = gamma * p_a / rho_a
        c_a = np.sqrt(c_a_sq)

        rho_s = self._znd.rho[-1]
        p_s = self._znd.p[-1]
        c_s = np.sqrt(self.c.gamma * p_s / rho_s)
        c_s_sq = c_s**2

        factors = {
            'density': rho_a / rho_s,
            'volume': rho_s / rho_a,
            'pressure': p_a / p_s,
            'velocity': c_a / (gamma_sqrt * c_s),
            'energy': c_a_sq / c_s_sq,
            'time': gamma_sqrt * c_s / c_a,
            'inv_time': c_a / (gamma_sqrt * c_s),
        }

        self._print_factors(factors, 'Lee & Stewart')

    def to_abouseif_toong_scales(self):
        k_AT = self._compute_k_in_abouseif_toong_scales()
        k_Erp = self._znd.k
        inv_t_factor = k_AT / k_Erp

        factors = {
            'time': 1.0 / inv_t_factor,
            'inv_time': inv_t_factor,
        }

        self._print_factors(factors, 'Abouseif & Toong')

    def _compute_k_in_abouseif_toong_scales(self):
        time, _ = integrate.quad(self._integrand_k_in_abouseif_toong_scales,
                                 0.0, 0.5)

        assert time > 0.0, 'Half-reaction time must be positive'

        return time

    def _integrand_k_in_abouseif_toong_scales(self, lamda):
        e_act, Q, g = self.c.e_act, self.c.q, self.c.gamma
        gp1 = g + 1

        M, P, H, h = self._znd.M, self._znd.P, self._znd.H, self._znd.h
        delta = np.sqrt(1 - h * M**2 * (H + Q * lamda) / P**2)
        v = g * P * (1 - delta) / (gp1 * M**2)
        p = P - M**2 * v

        numer = 1.0
        denom = (1 - lamda) * np.exp(-e_act / (p * v))
        result = numer / denom

        return result

    def _print_factors(self, factors, header):
        if self.trace:
            print('****************')
            print(header)

            for k, v in sorted(factors.items()):
                print('{:12s}: {:{fmt}}'.format(k, v, fmt=FORMAT))
