#!/usr/bin/env python
# -*- coding: utf-8 -*-

# Built-in imports
import itertools

# 3rd party imports
import numpy as np

__author__ = "Louis Richard"
__email__ = "louisr@irfu.se"
__copyright__ = "Copyright 2020-2023"
__license__ = "MIT"
__version__ = "2.4.2"
__status__ = "Prototype"


def _calc_diel(kc_, w_final, theta_, wp_e, wp_i, wc_i):
    # The elements of the dielectric tensor, using Swansons notation
    diel_s = 1 - wp_e**2 / (w_final**2 - 1) - wp_i**2 / (w_final**2 - wc_i**2)
    diel_d = -(wp_e**2) / (w_final * (w_final**2 - 1))
    diel_d += wc_i * wp_i**2 / (w_final * (w_final**2 - wc_i**2))
    diel_p = 1 - (wp_e**2 + wp_i**2) / w_final**2

    n2_ = kc_**2 / w_final**2

    diel_xx = diel_s - n2_ * np.cos(theta_) ** 2
    diel_xy = -1j * diel_d
    diel_xz = n2_ * np.cos(theta_) * np.sin(theta_)
    diel_yy = diel_s - n2_
    diel_zz = diel_p - n2_ * np.sin(theta_) ** 2
    return diel_xx, diel_xy, diel_xz, diel_yy, diel_zz


def _calc_e(diel_tensor):
    _, diel_xy, diel_xz, diel_yy, diel_zz = diel_tensor
    e_x = -diel_zz / diel_xz
    e_y = diel_xy / diel_yy * e_x
    e_z = np.ones(e_y.shape)

    e_per = np.sqrt(e_x * np.conj(e_x) + e_y * np.conj(e_y))
    e_pol = -2 * np.imag(e_x * np.conj(e_y)) / e_per**2
    e_tot = np.sqrt(e_x * np.conj(e_x) + e_y * np.conj(e_y) + e_z**2)

    return e_x, e_y, e_z, e_per, e_pol, e_tot


def _calc_b(kc_x_mat, kc_z_mat, w_final, e_x, e_y, e_z):
    b_x = -kc_z_mat * e_y / w_final
    b_y = (kc_z_mat * e_x - kc_x_mat * e_z) / w_final
    b_z = kc_x_mat * e_y / w_final

    b_par = np.sqrt(b_z * np.conj(b_z))
    b_per = np.sqrt(b_x * np.conj(b_x) + b_y * np.conj(b_y))
    b_pol = -2 * np.imag(b_x * np.conj(b_y)) / b_per**2
    b_tot = np.sqrt(
        b_x * np.conj(b_x) + b_y * np.conj(b_y) + b_z * np.conj(b_z),
    )

    return b_x, b_y, b_z, b_par, b_per, b_pol, b_tot


def _calc_s(e_x, e_y, e_z, b_x, b_y, b_z):
    # Poynting flux
    s_x = e_y * np.conj(b_z) - e_z * np.conj(b_y)
    s_y = e_z * np.conj(b_x) - e_x * np.conj(b_z)
    s_z = e_x * np.conj(b_y) - e_y * np.conj(b_x)
    s_par = np.abs(s_z)
    s_tot = np.sqrt(
        s_x * np.conj(s_x) + s_y * np.conj(s_y) + s_z * np.conj(s_z),
    )

    return s_par, s_tot


def _calc_part2fields(wp_e, en_e, en_i, e_tot, b_tot):
    n_e = wp_e**2
    en_e_n = en_e * n_e
    en_i_n = en_i * n_e
    en_efield = 0.5 * e_tot**2
    en_bfield = 0.5 * b_tot**2
    ratio_part_field = (en_e_n + en_i_n) / (en_efield + en_bfield)
    return ratio_part_field


def _calc_continuity(kc_x_mat, kc_z_mat, w_final, v_ex, v_ez, v_ix, v_iz):
    dn_e_n = (kc_x_mat * v_ex + kc_z_mat * v_ez) / w_final
    dn_e_n = np.sqrt(dn_e_n * np.conj(dn_e_n))
    dn_i_n = (kc_x_mat * v_ix + kc_z_mat * v_iz) / w_final
    dn_i_n = np.sqrt(dn_i_n * np.conj(dn_i_n))
    dne_dni = dn_e_n / dn_i_n

    return dn_e_n, dn_i_n, dne_dni


def _calc_vei(m_i, wc_i, w_final, e_x, e_y, e_z):
    q_e, q_i, m_e, wc_e = [-1, 1, 1, 1]

    v_ex = 1j * q_e * (w_final * e_x - 1j * wc_e * e_y)
    v_ex /= m_e * (w_final**2 - wc_e**2)

    v_ey = 1j * q_e * (1j * wc_e * e_x + w_final * e_y)
    v_ey /= m_e * (w_final**2 - wc_e**2)

    v_ez = 1j * q_e * e_z / (m_e * w_final)

    v_ix = 1j * q_i * (w_final * e_x + 1j * wc_i * e_y)
    v_ix /= m_i * (w_final**2 - wc_i**2)

    v_iy = 1j * q_i * (-1j * wc_i * e_x + w_final * e_y)
    v_iy /= m_i * (w_final**2 - wc_i**2)

    v_iz = 1j * q_i * e_z / (m_i * w_final)

    return v_ex, v_ey, v_ez, v_ix, v_iy, v_iz


def disp_surf_calc(kc_x_max, kc_z_max, m_i, wp_e):
    r"""Calculate the cold plasma dispersion surfaces according to equation
    2.64 in Plasma Waves by Swanson (2nd ed.)

    Parameters
    ----------
    kc_x_max : float
        Max value of k_perpendicular*c/w_c.
    kc_z_max : float
        Max value of k_parallel*c/w_c.
    m_i : float
        Ion mass in terms of electron masses.
    wp_e : float
        Electron plasma frequency in terms of electron gyro frequency.

    Returns
    -------
    kx_ : numpy.ndarray
        kperpandicular*c/w_c meshgrid
    kz_ : numpy.ndarray
        kparallel*c/w_c meshgrid
    wf_ : numpy.ndarray
        Dispersion surfaces.
    extra_param : dict
        Extra parameters to plot.
    """

    # Make vectors of the wave numbers
    kc_z = np.linspace(1e-6, kc_z_max, 35)
    kc_x = np.linspace(1e-6, kc_x_max, 35)

    # Turn those vectors into matrices
    kc_x_mat, kc_z_mat = np.meshgrid(kc_x, kc_z)

    # Find some of the numbers that appear later in the calculations
    kc_ = np.sqrt(kc_x_mat**2 + kc_z_mat**2)  # Absolute value of k
    theta_ = np.arctan2(kc_x_mat, kc_z_mat)  # The angle between k and B
    wc_i = 1 / m_i  # The ion gyro frequency
    wp_i = wp_e / np.sqrt(m_i)  # The ion plasma frequency
    wp_ = np.sqrt(wp_e**2 + wp_i**2)  # The total plasma frequency

    # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    # For every k_perp and k_par, turn the dispersion relation into a
    # polynomial equation and solve it.
    # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    # The polynomial coefficients are calculated
    pol_koeff_8 = -2 * kc_**2
    pol_koeff_8 -= (1 + wc_i**2 + 3 * wp_**2) * np.ones(kc_.shape)
    pol_koeff_6 = (2 * kc_**2 + wp_**2) * (1 + wc_i**2 + 2 * wp_**2)
    pol_koeff_6 += kc_**4 + (wp_**2 + wc_i) ** 2
    pol_koeff_4 = -(kc_**4) * (1 + wc_i**2 + wp_**2)
    pol_koeff_4 -= 2 * kc_**2 * (wp_**2 + wc_i) ** 2
    pol_koeff_4 -= (kc_ * wp_) ** 2 * (1 + wc_i**2 - wc_i) * (1 + np.cos(theta_) ** 2)
    pol_koeff_4 -= wp_**2 * (wp_**2 + wc_i) ** 2
    pol_koeff_2 = kc_**4 * (
        wp_**2 * (1 + wc_i**2 - wc_i) * np.cos(theta_) ** 2
        + wc_i * (wp_**2 + wc_i)
    )
    pol_koeff_2 += (
        kc_**2 * wp_**2 * wc_i * (wp_**2 + wc_i) * (1 + np.cos(theta_) ** 2)
    )
    pol_koeff_0 = -(kc_**4) * wc_i**2 * wp_**2 * np.cos(theta_) ** 2

    w_final = np.zeros((10, len(kc_z), len(kc_x)))

    # For each k, solve the equation
    for k_z, k_x in itertools.product(range(len(kc_z)), range(len(kc_x))):
        disp_polynomial = [
            1,
            0,
            pol_koeff_8[k_z, k_x],
            0,
            pol_koeff_6[k_z, k_x],
            0,
            pol_koeff_4[k_z, k_x],
            0,
            pol_koeff_2[k_z, k_x],
            0,
            pol_koeff_0[k_z, k_x],
        ]
        # theoretically should be real (A. Tjulin)
        w_temp = np.real(np.roots(disp_polynomial))
        # We need to sort the answers to get nice surfaces.
        w_final[:, k_z, k_x] = np.sort(w_temp)

    n2_ = kc_**2 / w_final**2
    v_ph_c = np.sqrt(1.0 / n2_)
    va_c = 1 / (wp_e * np.sqrt(m_i))
    v_ph_va = v_ph_c / va_c

    diel_tensor = _calc_diel(kc_, w_final, theta_, wp_e, wp_i, wc_i)

    e_x, e_y, e_z, _, e_pol, e_tot = _calc_e(diel_tensor)
    e_par = (kc_x_mat * e_x + kc_z_mat * e_z) / kc_

    b_x, b_y, b_z, b_par, _, b_pol, b_tot = _calc_b(
        kc_x_mat,
        kc_z_mat,
        w_final,
        e_x,
        e_y,
        e_z,
    )

    dk_x, dk_z = [kc_x_mat[1], kc_z_mat[1]]
    dw_x, dw_z = [np.zeros(w_final.shape) for _ in range(2)]
    dw_x[:, :, 1:] = np.diff(w_final, axis=2)
    dw_z[:, 1:, :] = np.diff(w_final, axis=1)
    v_x, v_z = [dw_ / dk for dw_, dk in zip([dw_x, dw_z], [dk_x, dk_z])]

    s_par, s_tot = _calc_s(e_x, e_y, e_z, b_x, b_y, b_z)

    # Compute ion and electron velocities
    v_ex, v_ey, v_ez, v_ix, v_iy, v_iz = _calc_vei(
        m_i,
        wc_i,
        w_final,
        e_x,
        e_y,
        e_z,
    )

    # Ratio of parallel and perpendicular to B speed
    vepar_perp = v_ez * np.conj(v_ez)
    vepar_perp /= v_ex * np.conj(v_ex) + v_ey * np.conj(v_ey)
    vipar_perp = v_iz * np.conj(v_iz)
    vipar_perp /= v_ix * np.conj(v_ix) + v_iy * np.conj(v_iy)

    # Total particle speeds
    v_e2 = v_ex * np.conj(v_ex) + v_ey * np.conj(v_ey) + v_ez * np.conj(v_ez)
    v_i2 = v_ix * np.conj(v_ix) + v_iy * np.conj(v_iy) + v_iz * np.conj(v_iz)

    # Ion and electron energies
    m_e = -1
    en_e = 0.5 * m_e * v_e2
    en_i = 0.5 * m_i * v_i2

    # Ratio of particle and field energy densities
    ratio_part_field = _calc_part2fields(wp_e, en_e, en_i, e_tot, b_tot)

    # Continuity equation
    dn_e_n, dn_i_n, dne_dni = _calc_continuity(
        kc_x_mat,
        kc_z_mat,
        w_final,
        v_ex,
        v_ez,
        v_ix,
        v_iz,
    )

    dn_e_n_db_b = dn_e_n / b_tot
    dn_i_n_db_b = dn_i_n / b_tot

    dn_e_n_dbpar_b = dn_e_n / b_par
    dn_i_n_dbpar_b = dn_i_n / b_par

    dn_e = dn_e_n * wp_e**2
    k_dot_e = e_x * kc_x_mat + e_z * kc_z_mat
    k_dot_e = np.sqrt(k_dot_e * np.conj(k_dot_e))

    # Build output dict
    extra_param = {
        "Degree of electromagnetism": np.log10(b_tot / e_tot),
        "Degree of longitudinality": np.abs(e_par) / e_tot,
        "Degree of parallelity E": e_z / e_tot,
        "Degree of parallelity B": np.sqrt(b_z * np.conj(b_z)) / b_tot,
        "W_E/W_B": np.log10(e_tot**2 / b_tot**2),
        "Ellipticity E": e_pol,
        "Ellipticity B": b_pol,
        "E_part/E_field": np.log10(ratio_part_field),
        "v_g": np.sqrt(v_x**2 + v_z**2),
        "v_ph/v_a": np.log10(v_ph_va),
        "E_e/E_i": np.log10(en_e / en_i),
        "v_e/v_i": np.log10(np.sqrt(v_e2 / v_i2)),
        "v_epara/v_eperp": np.log10(vepar_perp),
        "v_ipara/v_iperp": np.log10(vipar_perp),
        "dn_e/dn_i": np.log10(dne_dni),
        "(dn_e/n)/ (dB/B)": np.log10(dn_e_n_db_b),
        "(dn_i/n)/(dB/B)": np.log10(dn_i_n_db_b),
        "(dn_i/n)/(dBpar/B)": np.log10(dn_i_n_dbpar_b),
        "(dn_e/n)/(dB/B)": np.log10(dn_e / k_dot_e),
        "(dn_e/n)/(dBpar /B)": np.log10(dn_e_n_dbpar_b),
        " Spar/Stot": s_par / s_tot,
    }

    for k, v in zip(extra_param.keys(), extra_param.values()):
        extra_param[k] = np.transpose(np.real(v), [0, 2, 1])

    kx_ = np.transpose(kc_x_mat)
    kz_ = np.transpose(kc_z_mat)
    wf_ = np.transpose(w_final, [0, 2, 1])

    return kx_, kz_, wf_, extra_param
