"""
Integral computation for Bézier tariff curves.

Ported from eau.py:169-211 (NewModel.computeIntegrals).
Pure Python + numpy, no matplotlib.
"""

import numpy as np


def compute_integrals(
    volume: float,
    vinf: float,
    vmax: float,
    pmax: float,
    a: float,
    b: float,
    c: float,
    d: float,
    e: float,
) -> tuple[float, float, float]:
    """
    Compute (alpha1, alpha2, beta2) for a given consumption volume.

    The total bill for a household consuming `volume` m³ is:
        bill = abo + alpha1 * p0 + alpha2 * p0 + beta2

    where p0 is the inflection price (computed separately to balance revenue).

    Args:
        volume: consumption in m³ for this household
        vinf: inflection volume separating the two tiers
        vmax: maximum volume (price = pmax at this volume)
        pmax: maximum price per m³
        a, b: shape parameters for tier 1 Bézier curve
        c, d, e: shape parameters for tier 2 Bézier curve

    Returns:
        (alpha1, alpha2, beta2) tuple
    """
    if volume <= vinf:
        # Tier 1 only
        T = _solve_tier1_t(volume, vinf, b)
        alpha1 = _compute_alpha1(T, vinf, a, b)
        return alpha1, 0.0, 0.0
    else:
        # Full tier 1 (T=1) + partial tier 2
        alpha1 = _compute_alpha1(1.0, vinf, a, b)

        # Tier 2
        wmax = vmax - vinf
        T = _solve_tier2_t(volume - vinf, wmax, c, d)

        uu = _compute_uu(T, c, d, e)
        alpha2 = (volume - vinf) - 3 * uu * wmax
        beta2 = 3 * pmax * wmax * uu
        return alpha1, alpha2, beta2


def _solve_tier1_t(volume: float, vinf: float, b: float) -> float:
    """Find T such that v(T) = volume for tier 1."""
    if volume == 0:
        return 0.0
    if volume >= vinf:
        return 1.0

    # Solve: vinf * [(1 - 3b) * T³ + 3b * T²] = volume
    # => (1-3b) * T³ + 3b * T² - volume/vinf = 0
    p = [1 - 3 * b, 3 * b, 0, -volume / vinf]
    roots = np.roots(p)
    roots = np.unique(roots)
    real_roots = np.real(roots[np.isreal(roots)])
    mask = (real_roots <= 1.0) & (real_roots >= 0.0)
    return float(real_roots[mask][0])


def _solve_tier2_t(w: float, wmax: float, c: float, d: float) -> float:
    """Find T such that w(T) = w for tier 2, where w = volume - vinf."""
    if w == 0:
        return 0.0
    if w >= wmax:
        return 1.0

    # Solve: wmax * [(3(c+d-cd)-2)*T³ + 3(1-2c-d+cd)*T² + 3c*T] = w
    p = [
        3 * (c + d - c * d) - 2,
        3 * (1 - 2 * c - d + c * d),
        3 * c,
        -w / wmax,
    ]
    roots = np.roots(p)
    roots = np.unique(roots)
    real_roots = np.real(roots[np.isreal(roots)])
    mask = (real_roots <= 1.0 + 1e-10) & (real_roots >= -1e-10)
    if not mask.any():
        # Fallback: closest root to [0,1]
        return float(np.clip(np.real(roots[0]), 0.0, 1.0))
    return float(np.clip(real_roots[mask][0], 0.0, 1.0))


def _compute_alpha1(T: float, vinf: float, a: float, b: float) -> float:
    """Compute alpha1 coefficient for tier 1."""
    return 3 * vinf * (
        T**6 / 6 * (-9 * a * b + 3 * a + 6 * b - 2)
        + T**5 / 5 * (24 * a * b - 6 * a - 13 * b + 3)
        + 3 * T**4 / 4 * (-7 * a * b + a + 2 * b)
        + T**3 / 3 * 6 * a * b
    )


def _compute_uu(T: float, c: float, d: float, e: float) -> float:
    """Compute the uu intermediate value for tier 2."""
    return (
        (-3 * c * d + 9 * e * c * d + 3 * c - 9 * e * c + 3 * d - 9 * e * d + 6 * e - 2) * T**6 / 6
        + (2 * c * d - 15 * e * c * d - 4 * c + 21 * e * c - 2 * d + 15 * e * d - 12 * e + 2) * T**5 / 5
        + (6 * e * c * d + c - 15 * e * c - 6 * e * d + 6 * e) * T**4 / 4
        + (3 * e * c) * T**3 / 3
    )