"""
Shannon Equation - Core Calculations Module

All scientific computation functions extracted from the original Shannon.py.
These are pure functions with no UI dependency.
"""

from math import log, pi, sqrt, cos, acos, atan
import numpy as np

# ──────────────────────────────────────────────────────────────────────────────
#  Fundamental Shannon Functions
# ──────────────────────────────────────────────────────────────────────────────

def combine_cnr(*cnr_values: float) -> float:
    """Combine multiple Carrier-to-Noise Ratios (in dB) into one equivalent C/N.

    Uses the summation of normalized noise variances:
        1/CNR_total = sum(1/CNR_i)
    """
    ncr_linear = 0.0
    for cnr_db in cnr_values:
        ncr_linear += 10 ** (-cnr_db / 10)
    return -10 * log(ncr_linear, 10)


def shannon_capacity(bw: float = 36.0, cnr: float = 10.0, penalty: float = 0.0) -> float:
    """Shannon channel capacity (bit rate in Mbps).

    Args:
        bw: Bandwidth in MHz.
        cnr: Carrier-to-Noise Ratio in dB.
        penalty: Implementation penalty in dB.
    """
    cnr_linear = 10 ** ((cnr - penalty) / 10)
    return bw * log(1 + cnr_linear, 2)


def br_multiplier(bw_mul: float = 1.0, p_mul: float = 2.0, cnr: float = 10.0) -> float:
    """Bit Rate multiplying factor when BW and Power are scaled."""
    cnr_linear = 10 ** (cnr / 10)
    return bw_mul * log(1 + cnr_linear * p_mul / bw_mul, 2) / log(1 + cnr_linear, 2)


def shannon_points(bw: float = 36.0, cnr: float = 10.0):
    """Compute key Shannon operating points.

    Returns:
        (cnr_linear, br_infinity, c_n0, br_constrained)
    """
    cnr_linear = 10 ** (cnr / 10)
    c_n0 = cnr_linear * bw
    br_infinity = c_n0 / log(2)
    br_constrained = shannon_capacity(bw, cnr)
    return cnr_linear, br_infinity, c_n0, br_constrained


# ──────────────────────────────────────────────────────────────────────────────
#  Satellite Link Budget Calculations
# ──────────────────────────────────────────────────────────────────────────────

def compute_satellite_link(
    freq_ghz: float,
    hpa_power_w: float,
    sat_loss_db: float,
    sat_cir_list: list[float],
    sat_beam_deg: float,
    gain_offset_db: float,
    sat_alt_km: float,
    sat_lat: float,
    sat_lon: float,
    gs_lat: float,
    gs_lon: float,
    availability_pct: float,
) -> dict:
    """Compute all satellite link parameters.

    Returns a dict with all computed values.
    """
    import itur

    R_EARTH = 6378  # km
    SAT_ANT_EFF = 0.65

    # Signal power after losses
    sig_power = hpa_power_w * 10 ** (-sat_loss_db / 10)

    # Satellite antenna
    wavelength = 300e6 / freq_ghz / 1e9  # meters
    sat_gain_linear = SAT_ANT_EFF * (pi * 70 / sat_beam_deg) ** 2
    sat_gain_linear *= 10 ** (-gain_offset_db / 10)

    # EIRP
    eirp_linear = sig_power * sat_gain_linear

    # Satellite C/I
    sat_cir = combine_cnr(*sat_cir_list)

    # Path geometry
    path_length = sqrt(
        sat_alt_km ** 2
        + 2 * R_EARTH * (R_EARTH + sat_alt_km)
        * (1 - cos(np.radians(sat_lat - gs_lat)) * cos(np.radians(sat_lon - gs_lon)))
    )

    phi = acos(cos(np.radians(sat_lat - gs_lat)) * cos(np.radians(sat_lon - gs_lon)))
    if phi > 0:
        elevation = float(np.degrees(
            atan((cos(phi) - R_EARTH / (R_EARTH + sat_alt_km)) / sqrt(1 - cos(phi) ** 2))
        ))
    else:
        elevation = 90.0

    # Atmospheric attenuation
    if elevation <= 0:
        atm_loss_db = 999.0
    else:
        atm_loss_db = float(
            itur.atmospheric_attenuation_slant_path(
                gs_lat, gs_lon, freq_ghz, elevation, 100 - availability_pct, 1
            ).value
        )

    # Path dispersion
    path_loss_linear = 4 * pi * (path_length * 1000) ** 2
    free_space_loss_linear = (4 * pi * path_length * 1000 / wavelength) ** 2

    # PFD
    pfd_linear = eirp_linear / path_loss_linear * 10 ** (-atm_loss_db / 10)

    return {
        "sig_power": sig_power,
        "wavelength": wavelength,
        "sat_gain_linear": sat_gain_linear,
        "eirp_linear": eirp_linear,
        "sat_cir": sat_cir,
        "path_length": path_length,
        "elevation": elevation,
        "atm_loss_db": atm_loss_db,
        "path_loss_linear": path_loss_linear,
        "free_space_loss_linear": free_space_loss_linear,
        "pfd_linear": pfd_linear,
    }


def compute_receiver(
    pfd_linear: float,
    atm_loss_db: float,
    wavelength: float,
    cpe_ant_d: float,
    cpe_t_clear: float,
) -> dict:
    """Compute receiver-side parameters."""
    CPE_ANT_EFF = 0.6
    K_BOLTZ = 1.38e-23  # J/K

    cpe_t_att = (cpe_t_clear - 40) + 40 * 10 ** (-atm_loss_db / 10) + 290 * (1 - 10 ** (-atm_loss_db / 10))

    cpe_ae = pi * cpe_ant_d ** 2 / 4 * CPE_ANT_EFF
    cpe_gain_linear = (pi * cpe_ant_d / wavelength) ** 2 * CPE_ANT_EFF
    cpe_g_t = 10 * log(cpe_gain_linear / cpe_t_att, 10)

    rx_power = pfd_linear * cpe_ae
    n0 = K_BOLTZ * cpe_t_att
    c_n0_hz = rx_power / n0
    c_n0_mhz = c_n0_hz / 1e6

    br_infinity = c_n0_mhz / log(2)

    # Spectral efficiency points
    bw_spe_1 = c_n0_mhz
    bw_spe_double = c_n0_mhz / (2 ** 2 - 1)

    br_spe_1 = bw_spe_1
    br_spe_double = bw_spe_double * 2

    return {
        "cpe_ae": cpe_ae,
        "cpe_gain_linear": cpe_gain_linear,
        "cpe_g_t": cpe_g_t,
        "cpe_t_att": cpe_t_att,
        "rx_power": rx_power,
        "n0": n0,
        "c_n0_hz": c_n0_hz,
        "c_n0_mhz": c_n0_mhz,
        "br_infinity": br_infinity,
        "bw_spe_1": bw_spe_1,
        "br_spe_1": br_spe_1,
        "br_spe_double": br_spe_double,
    }


def compute_baseband(
    c_n0_mhz: float,
    br_infinity: float,
    bw_spe_1: float,
    sat_cir: float,
    bandwidth: float,
    rolloff: float,
    overheads: float,
    cnr_imp_list: list[float],
    penalties: float,
) -> dict:
    """Compute baseband processing results."""
    cnr_imp = combine_cnr(*cnr_imp_list)

    cnr_spe_1 = 0.0  # dB
    cnr_bw = cnr_spe_1 + 10 * log(bw_spe_1 / bandwidth, 10)
    bw_nyq = bandwidth / (1 + rolloff / 100)
    cnr_nyq = cnr_spe_1 + 10 * log(bw_spe_1 / bw_nyq, 10)
    cnr_rcv = combine_cnr(cnr_nyq, cnr_imp, sat_cir)

    br_nyq = shannon_capacity(bw_nyq, cnr_nyq)
    br_rcv = shannon_capacity(bw_nyq, cnr_rcv, penalties)
    br_rcv_higher = br_rcv / (1 + overheads / 100)

    spe_nyq = br_nyq / bandwidth
    bits_per_symbol = br_nyq / bw_nyq
    spe_rcv = br_rcv / bandwidth
    spe_higher = br_rcv_higher / bandwidth

    return {
        "cnr_bw": cnr_bw,
        "cnr_nyq": cnr_nyq,
        "cnr_rcv": cnr_rcv,
        "cnr_imp": cnr_imp,
        "bw_nyq": bw_nyq,
        "br_nyq": br_nyq,
        "br_rcv": br_rcv,
        "br_rcv_higher": br_rcv_higher,
        "br_nyq_norm": br_nyq / br_infinity,
        "br_rcv_norm": br_rcv / br_infinity,
        "br_rcv_h_norm": br_rcv_higher / br_infinity,
        "spe_nyq": spe_nyq,
        "bits_per_symbol": bits_per_symbol,
        "spe_rcv": spe_rcv,
        "spe_higher": spe_higher,
        "bandwidth": bandwidth,
        "rolloff": rolloff,
        "overheads": overheads,
        "penalties": penalties,
    }


# ──────────────────────────────────────────────────────────────────────────────
#  Formatting Helpers
# ──────────────────────────────────────────────────────────────────────────────

def fmt_br(br: float) -> str:
    return f"{br:.1f} Mbps"


def fmt_power(p: float) -> str:
    p_db = 10 * log(p, 10)
    if 1 < p < 1e4:
        return f"{p:.1f} W · {p_db:.1f} dBW"
    elif 1e-3 < p <= 1:
        return f"{p:.4f} W · {p_db:.1f} dBW"
    else:
        return f"{p:.1e} W · {p_db:.1f} dBW"


def fmt_pfd(p: float) -> str:
    p_db = 10 * log(p, 10)
    return f"{p:.1e} W/m² · {p_db:.1f} dBW/m²"


def fmt_psd(p: float) -> str:
    p_db = 10 * log(p, 10)
    return f"{p:.1e} W/MHz · {p_db:.1f} dBW/MHz"


def fmt_gain(g: float) -> str:
    g_db = 10 * log(g, 10)
    return f"{g:.1f} · {g_db:.1f} dBi"


def fmt_ploss(loss: float) -> str:
    loss_db = 10 * log(loss, 10)
    return f"{loss:.2e} m² · {loss_db:.1f} dBm²"