#  Copyright (c) 2026. Jose M. Requena-Plens
"""Tests for ISO/TS 7849 airborne sound power from surface vibration.

Anchored on the standard's own worked calibration example (Part 1, Eq. 8:
a_peak = 9,81 m/s^2 at 100 Hz -> L_v = 106,9 dB), on the closed-form velocity
and sound-power-level relations (Eqs. 3, 12/15) and on the exact round-trip
between the radiation factor (Eq. 8: eps = P/(Z_c v^2 S)) and the sound power
level L_W = 10 lg(P/P0). The K1A extraneous-velocity correction reproduces
Table 2 exactly.
"""

from __future__ import annotations

import math

import numpy as np
import pytest

from phonometry import (
    VibrationSoundPowerResult,
    extraneous_velocity_correction,
    mean_velocity_level,
    radiated_sound_power_level,
    radiation_factor,
    sound_power_from_vibration,
    velocity_level,
    velocity_level_from_acceleration,
)

V0 = 5.0e-8
ZCN, ZC0, P0 = 411.0, 400.0, 1.0e-12


# ---------------------------------------------------------------------------
# Velocity level (Eq. 3) and calibration (Eq. 8)
# ---------------------------------------------------------------------------
def test_velocity_level_reference() -> None:
    assert velocity_level(V0) == pytest.approx(0.0)
    assert velocity_level(5.0e-5) == pytest.approx(60.0)   # 20 lg(1000)


def test_calibration_example_from_standard() -> None:
    """ISO/TS 7849-1 worked EXAMPLE: 9,81 m/s^2 at 100 Hz -> 106,9 dB."""
    lv = float(velocity_level_from_acceleration(9.81, 100.0))
    assert lv == pytest.approx(106.9, abs=0.05)


def test_calibration_matches_hand_formula() -> None:
    a, f = 4.0, 250.0
    expected = 20.0 * math.log10(a / (2.0 * math.pi * f * V0 * math.sqrt(2.0)))
    assert velocity_level_from_acceleration(a, f) == pytest.approx(expected)


def test_calibration_rejects_non_positive_frequency() -> None:
    with pytest.raises(ValueError, match="frequency"):
        velocity_level_from_acceleration(9.81, 0.0)


# ---------------------------------------------------------------------------
# Mean velocity level (Eq. 10/11)
# ---------------------------------------------------------------------------
def test_mean_velocity_level_energetic() -> None:
    levels = np.array([60.0, 60.0, 66.0])
    expected = 10.0 * math.log10(np.mean(10.0 ** (0.1 * levels)))
    assert mean_velocity_level(levels) == pytest.approx(expected)


def test_mean_velocity_level_area_weighted() -> None:
    levels = np.array([60.0, 66.0])
    areas = np.array([1.0, 3.0])
    expected = 10.0 * math.log10(
        np.sum(areas * 10.0 ** (0.1 * levels)) / np.sum(areas)
    )
    assert mean_velocity_level(levels, areas=areas) == pytest.approx(expected)


def test_mean_velocity_level_area_shape_mismatch() -> None:
    with pytest.raises(ValueError, match="areas"):
        mean_velocity_level([60.0, 66.0], areas=[1.0])


# ---------------------------------------------------------------------------
# Radiation factor (Eq. 4/8) and the round-trip to the sound power level
# ---------------------------------------------------------------------------
def test_radiation_factor_definition() -> None:
    p, s, v2 = 3.0e-4, 2.0, (1.0e-3) ** 2
    assert radiation_factor(p, s, v2) == pytest.approx(p / (ZCN * v2 * s))


def test_power_level_round_trip_through_radiation_factor() -> None:
    """eps from Eq. 8, fed into Eq. 15, recovers L_W = 10 lg(P/P0) exactly."""
    p, s, v2 = 3.0e-4, 2.0, (1.0e-3) ** 2
    eps = float(radiation_factor(p, s, v2))
    lv = float(velocity_level(math.sqrt(v2)))
    lw = float(radiated_sound_power_level(lv, s, radiation_factor=eps))
    assert lw == pytest.approx(10.0 * math.log10(p / P0))


def test_power_level_impedance_constant() -> None:
    # eps = 1, S = S0 = 1 -> L_W = L_v + 10 lg(411/400)
    lv = 80.0
    lw = float(radiated_sound_power_level(lv, 1.0))
    assert lw == pytest.approx(lv + 10.0 * math.log10(ZCN / ZC0))


def test_upper_limit_is_largest() -> None:
    # Part 1 (eps = 1) is an upper limit for any eps < 1.
    lv, s = 75.0, 2.0
    upper = float(radiated_sound_power_level(lv, s))
    measured = float(radiated_sound_power_level(lv, s, radiation_factor=0.4))
    assert upper > measured


def test_power_level_rejects_bad_area() -> None:
    with pytest.raises(ValueError, match="area"):
        radiated_sound_power_level(80.0, 0.0)


# ---------------------------------------------------------------------------
# Extraneous-velocity correction (Table 2)
# ---------------------------------------------------------------------------
def test_k1a_table_values() -> None:
    table = {3: 3.0, 4: 2.0, 5: 2.0, 6: 1.0, 7: 1.0, 8: 1.0, 9: 1.0, 10: 0.0}
    for dlv, k in table.items():
        assert extraneous_velocity_correction(float(dlv)) == pytest.approx(k)


def test_k1a_boundaries() -> None:
    assert extraneous_velocity_correction(2.0) == 3.0    # dLv < 3 -> 3 dB
    assert extraneous_velocity_correction(15.0) == 0.0   # dLv >= 10 -> 0


# ---------------------------------------------------------------------------
# Result bundle
# ---------------------------------------------------------------------------
def test_result_bundle_and_total() -> None:
    lv = np.array([70.0, 75.0, 72.0])
    eps = np.array([0.5, 0.8, 1.0])
    f = np.array([250.0, 500.0, 1000.0])
    res = sound_power_from_vibration(lv, 1.5, radiation_factor=eps, frequencies=f)
    assert isinstance(res, VibrationSoundPowerResult)
    assert res.area == 1.5
    # per-band level matches the standalone function
    assert np.allclose(
        res.sound_power_level,
        radiated_sound_power_level(lv, 1.5, radiation_factor=eps),
    )
    # total is the energetic sum of the bands
    expected_total = 10.0 * math.log10(np.sum(10.0 ** (0.1 * res.sound_power_level)))
    assert res.total_level == pytest.approx(expected_total)


def test_result_frequencies_shape_mismatch() -> None:
    with pytest.raises(ValueError, match="frequencies"):
        sound_power_from_vibration([70.0, 75.0, 72.0], 1.5, frequencies=[250.0, 500.0])


def test_result_scalar_frequency_is_coerced() -> None:
    # a single-band call with a scalar frequency plots without error
    res = sound_power_from_vibration(80.0, 2.0, frequencies=1000.0)
    assert res.frequencies is not None
    assert res.frequencies.shape == (1,)


def test_result_scalar_radiation_factor_broadcasts() -> None:
    res = sound_power_from_vibration(np.array([70.0, 72.0]), 2.0)
    assert res.radiation_factor.shape == (2,)
    assert np.all(res.radiation_factor == 1.0)


def test_plot_returns_axes() -> None:
    pytest.importorskip("matplotlib")
    import matplotlib

    matplotlib.use("Agg")
    res = sound_power_from_vibration(
        np.array([70.0, 75.0, 72.0]), 1.5, frequencies=np.array([250.0, 500.0, 1000.0])
    )
    assert res.plot() is not None
