#  Copyright (c) 2026. Jose M. Requena-Plens
"""
Normal equal-loudness-level contours per ISO 226:2023.

Implements Formula (1) (clause 4.1: SPL of a pure tone from its loudness
level) and Formula (2) (clause 4.2: the inverse) with the Table 1 (p. 4)
parameters at the 29 preferred third-octave frequencies of ISO 266.
"""

from __future__ import annotations

from typing import Tuple

import numpy as np

# ISO 226:2023 Table 1 (p. 4): frequency Hz -> (alpha_f, L_U dB, T_f dB).
# alpha_f: exponent for loudness perception; L_U: magnitude of the linear
# transfer function normalized at 1 kHz; T_f: threshold of hearing.
_TABLE1: dict[float, tuple[float, float, float]] = {
    20.0: (0.635, -31.5, 78.1),
    25.0: (0.602, -27.2, 68.7),
    31.5: (0.569, -23.1, 59.5),
    40.0: (0.537, -19.3, 51.1),
    50.0: (0.509, -16.1, 44.0),
    63.0: (0.482, -13.1, 37.5),
    80.0: (0.456, -10.4, 31.5),
    100.0: (0.433, -8.2, 26.5),
    125.0: (0.412, -6.3, 22.1),
    160.0: (0.391, -4.6, 17.9),
    200.0: (0.373, -3.2, 14.4),
    250.0: (0.357, -2.1, 11.4),
    315.0: (0.343, -1.2, 8.6),
    400.0: (0.330, -0.5, 6.2),
    500.0: (0.320, 0.0, 4.4),
    630.0: (0.311, 0.4, 3.0),
    800.0: (0.303, 0.5, 2.2),
    1000.0: (0.300, 0.0, 2.4),
    1250.0: (0.295, -2.7, 3.5),
    1600.0: (0.292, -4.2, 1.7),
    2000.0: (0.290, -1.2, -1.3),
    2500.0: (0.290, 1.4, -4.2),
    3150.0: (0.289, 2.3, -6.0),
    4000.0: (0.289, 1.0, -5.4),
    5000.0: (0.289, -2.3, -1.5),
    6300.0: (0.293, -7.2, 6.0),
    8000.0: (0.303, -11.2, 12.6),
    10000.0: (0.323, -10.9, 13.9),
    12500.0: (0.354, -3.5, 12.3),
}

_FREQUENCIES = np.array(sorted(_TABLE1))


def _params(frequency: float) -> tuple[float, float, float]:
    """Table 1 parameters for a preferred third-octave frequency.

    The 1e-6 relative tolerance only absorbs floating-point representation
    error; adjacent table rows are >20 % apart, so a near-miss can never
    resolve to the wrong row - anything else raises.
    """
    for f, params in _TABLE1.items():
        if abs(frequency - f) <= 1e-6 * f:
            return params
    raise ValueError(
        f"ISO 226:2023 Table 1 defines parameters only at the 29 preferred "
        f"third-octave frequencies (20 Hz to 12.5 kHz); got frequency "
        f"{frequency!r} Hz. The standard specifies no interpolation."
    )


def _spl_from_phon(frequency: float, phon: float) -> float:
    """ISO 226:2023 Formula (1), clause 4.1 (p. 2)."""
    alpha_f, l_u, t_f = _params(frequency)
    term = (4.0e-10) ** (0.3 - alpha_f) * (
        10 ** (0.03 * phon) - 10 ** 0.072
    ) + 10 ** (alpha_f * (t_f + l_u) / 10)
    return float(10 / alpha_f * np.log10(term) - l_u)


def equal_loudness_contour(phon: float) -> Tuple[np.ndarray, np.ndarray]:
    """
    Normal equal-loudness-level contour (ISO 226:2023 Formula 1).

    Returns the sound pressure levels of pure tones judged equally loud as
    a 1 kHz tone at ``phon`` dB SPL, at the 29 preferred third-octave
    frequencies of Table 1.

    Validity (clause 4.1): 20 phon to 90 phon between 20 Hz and 4 kHz, and
    up to 80 phon between 5 kHz and 12.5 kHz — above 80 phon the returned
    contour therefore stops at 4 kHz.

    :param phon: Loudness level in phons (20 to 90).
    :return: Tuple ``(frequencies, spl)`` in Hz and dB re 20 uPa.
    """
    if not 20.0 <= phon <= 90.0:
        raise ValueError(
            "ISO 226:2023 Formula (1) is specified for 20 phon to 90 phon "
            f"(80 phon above 4 kHz); got {phon!r}."
        )
    freqs = _FREQUENCIES if phon <= 80.0 else _FREQUENCIES[_FREQUENCIES <= 4000.0]
    spl = np.array([_spl_from_phon(f, phon) for f in freqs])
    return freqs.copy(), spl


def loudness_level(spl: float, frequency: float) -> float:
    """
    Loudness level of a pure tone (ISO 226:2023 Formula 2).

    :param spl: Sound pressure level of the tone in dB re 20 uPa.
    :param frequency: Tone frequency in Hz; must be one of the 29 preferred
        third-octave frequencies of Table 1 (the standard specifies no
        interpolation between them).
    :return: Loudness level in phons. Values outside 20-90 phon (80 above
        4 kHz) are extrapolations the standard labels as informative only.
    """
    alpha_f, l_u, t_f = _params(frequency)
    b = (
        10 ** (alpha_f * (spl + l_u) / 10) - 10 ** (alpha_f * (t_f + l_u) / 10)
    ) / (4.0e-10) ** (0.3 - alpha_f) + 10 ** 0.072
    return float(100.0 / 3.0 * np.log10(b))


def hearing_threshold() -> Tuple[np.ndarray, np.ndarray]:
    """
    Threshold of hearing T_f (ISO 226:2023 Table 1).

    :return: Tuple ``(frequencies, threshold)`` in Hz and dB re 20 uPa.
    """
    tf = np.array([_TABLE1[f][2] for f in _FREQUENCIES])
    return _FREQUENCIES.copy(), tf
