#  Copyright (c) 2026. Jose M. Requena-Plens
"""
Installed structure-borne sound from service equipment (EN 12354-5:2009).

EN 12354-5 predicts the sound pressure level in a receiving room caused by
building service equipment that injects **structure-borne sound** into the
building. The chain closes the structural-vibroacoustics series:

1. The source strength is its *characteristic structure-borne sound power level*
   ``L_Ws,c`` (measured per EN 15657; see :mod:`phonometry.structure_borne_power`).
2. Only part of that power is actually injected into the supporting element; the
   loss is the **coupling term** ``D_C`` (clause 4.4.3), always positive, set by
   the source mobility ``Y_s`` and the receiver mobility ``Y_i`` (Formula 19b)::

       D_C,i = 10 lg( |Y_s + Y_i|**2 / (|Y_s| Re{Y_i}) )

   which reduces to ``10 lg(|Y_s|/Re{Y_i})`` for a force source (high source
   mobility, Formula 19c) and to ``-10 lg(|Y_s| Re{Z_i})`` for a velocity source
   (low source mobility, Formula 19d). An elastic support adds its transfer
   mobility ``Y_k`` inside the modulus (Formula 19e).
3. The **installed** power level is then (Formula 18b)::

       L_Ws,inst,i = L_Ws,c - D_C,i

4. The normalised sound pressure level in the receiving room for one path (i->j)
   follows from the installed power, the structure-to-airborne adjustment term
   ``D_sa`` (clause 4.4.4), the flanking sound reduction index ``R_ij,ref`` and
   the element area (Formula 18a)::

       L_n,s,ij = L_Ws,inst,i - D_sa,i - R_ij,ref - 10 lg(S_i/S0) - 10 lg(A0/4)

   with ``S0 = A0 = 10 m2``; the paths combine energetically (Formula 17).

The source and receiver mobilities/impedances are those of
:mod:`phonometry.mechanical_mobility` and :mod:`phonometry.transfer_stiffness`.
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import TYPE_CHECKING, Any

import numpy as np

if TYPE_CHECKING:
    from matplotlib.axes import Axes

from numpy.typing import ArrayLike

from ._validation import require_positive

#: Reference area ``S0 = A0`` of EN 12354-5 (Formula 18a), m^2.
REFERENCE_AREA: float = 10.0


def coupling_term(
    source_mobility: ArrayLike,
    receiver_mobility: ArrayLike,
    *,
    transfer_mobility: ArrayLike = 0.0,
) -> np.ndarray:
    """Coupling term ``D_C`` for a point excitation (EN 12354-5, Formula 19b/19e).

    ``D_C = 10 lg( |Y_s + Y_i + Y_k|**2 / (|Y_s| Re{Y_i}) )`` -- the loss between
    the characteristic and the injected structure-borne power. ``Y_k`` is the
    transfer mobility of an elastic support (Formula 19e; 0 for a rigid
    connection, Formula 19b).

    :param source_mobility: Source point mobility ``Y_s`` (complex), in m/(N.s).
    :param receiver_mobility: Receiver point mobility ``Y_i`` (complex).
    :param transfer_mobility: Elastic-support transfer mobility ``Y_k``
        (Default: 0.0).
    :return: The coupling term ``D_C``, in dB (>= 0 for passive systems).
    """
    ys = np.asarray(source_mobility, dtype=np.complex128)
    yi = np.asarray(receiver_mobility, dtype=np.complex128)
    yk = np.asarray(transfer_mobility, dtype=np.complex128)
    numerator = np.abs(ys + yi + yk) ** 2
    denominator = np.abs(ys) * np.real(yi)
    return np.asarray(10.0 * np.log10(numerator / denominator), dtype=np.float64)


def coupling_term_force_source(
    source_mobility: ArrayLike, receiver_mobility: ArrayLike
) -> np.ndarray:
    """Coupling term for a force source, high source mobility (Formula 19c).

    ``D_C = 10 lg(|Y_s| / Re{Y_i})``.

    :param source_mobility: Source point mobility ``Y_s`` (complex).
    :param receiver_mobility: Receiver point mobility ``Y_i`` (complex).
    :return: The coupling term ``D_C``, in dB.
    """
    ys = np.abs(np.asarray(source_mobility, dtype=np.complex128))
    yi = np.real(np.asarray(receiver_mobility, dtype=np.complex128))
    return np.asarray(10.0 * np.log10(ys / yi), dtype=np.float64)


def coupling_term_velocity_source(
    source_mobility: ArrayLike, receiver_impedance: ArrayLike
) -> np.ndarray:
    """Coupling term for a velocity source, low source mobility (Formula 19d).

    ``D_C = -10 lg(|Y_s| Re{Z_i})``.

    :param source_mobility: Source point mobility ``Y_s`` (complex).
    :param receiver_impedance: Receiver point impedance ``Z_i`` (complex).
    :return: The coupling term ``D_C``, in dB.
    """
    ys = np.abs(np.asarray(source_mobility, dtype=np.complex128))
    zi = np.real(np.asarray(receiver_impedance, dtype=np.complex128))
    return np.asarray(-10.0 * np.log10(ys * zi), dtype=np.float64)


def installed_structure_borne_power_level(
    characteristic_power_level: ArrayLike, coupling_term: ArrayLike
) -> np.ndarray:
    """Installed structure-borne power level (EN 12354-5, Formula 18b).

    ``L_Ws,inst,i = L_Ws,c - D_C,i``.

    :param characteristic_power_level: Characteristic level ``L_Ws,c`` (per
        band), in dB (from EN 15657).
    :param coupling_term: Coupling term ``D_C,i`` (per band), in dB.
    :return: The installed structure-borne power level ``L_Ws,inst``, in dB.
    """
    lw = np.asarray(characteristic_power_level, dtype=np.float64)
    dc = np.asarray(coupling_term, dtype=np.float64)
    return np.asarray(lw - dc, dtype=np.float64)


def structure_borne_pressure_level_path(
    installed_power_level: ArrayLike,
    adjustment_term: ArrayLike,
    flanking_reduction_index: ArrayLike,
    element_area: float,
    *,
    reference_area: float = REFERENCE_AREA,
) -> np.ndarray:
    """Normalised structure-borne SPL for one path i->j (Formula 18a).

    ``L_n,s,ij = L_Ws,inst,i - D_sa,i - R_ij,ref - 10 lg(S_i/S0) - 10 lg(A0/4)``.

    :param installed_power_level: Installed power level ``L_Ws,inst,i``, in dB.
    :param adjustment_term: Structure-to-airborne adjustment ``D_sa,i`` (clause
        4.4.4 / Annex F), in dB.
    :param flanking_reduction_index: Flanking sound reduction index
        ``R_ij,ref`` re ``S0`` (EN 12354-1), in dB.
    :param element_area: Supporting-element area ``S_i``, in m^2 (> 0).
    :param reference_area: Reference area ``S0 = A0`` (Default: 10 m^2).
    :return: The normalised path sound pressure level ``L_n,s,ij``, in dB.
    :raises ValueError: for a non-positive area.
    """
    element_area = require_positive(element_area, "element_area")
    reference_area = require_positive(reference_area, "reference_area")
    lw = np.asarray(installed_power_level, dtype=np.float64)
    dsa = np.asarray(adjustment_term, dtype=np.float64)
    rij = np.asarray(flanking_reduction_index, dtype=np.float64)
    lp = (
        lw - dsa - rij
        - 10.0 * np.log10(element_area / reference_area)
        - 10.0 * np.log10(reference_area / 4.0)
    )
    return np.asarray(lp, dtype=np.float64)


def total_structure_borne_pressure_level(path_levels: ArrayLike) -> np.ndarray:
    """Combine path sound pressure levels energetically (Formula 17).

    ``L_n,s = 10 lg( sum_j 10^(L_n,s,ij/10) )``.

    :param path_levels: Path levels ``L_n,s,ij``; sum is over the first axis
        (paths), broadcasting any trailing band axis.
    :return: The total normalised sound pressure level ``L_n,s``, in dB.
    """
    lp = np.asarray(path_levels, dtype=np.float64)
    return np.asarray(
        10.0 * np.log10(np.sum(10.0 ** (0.1 * lp), axis=0)), dtype=np.float64
    )


@dataclass(frozen=True)
class InstalledSourceResult:
    """Installed structure-borne sound prediction (EN 12354-5).

    :ivar frequencies: Band centre frequencies, in hertz, or ``None``.
    :ivar path_levels: Per-path normalised SPL ``L_n,s,ij`` (paths x bands), dB.
    :ivar total_level: Combined normalised SPL ``L_n,s`` per band, in dB.
    :ivar installed_power_level: Installed power level ``L_Ws,inst`` per band, dB.
    """

    path_levels: np.ndarray
    total_level: np.ndarray
    installed_power_level: np.ndarray
    frequencies: np.ndarray | None = None

    @property
    def overall_level(self) -> float:
        """Band-summed total level ``10 lg(sum 10^(0.1 L_n,s))``, in dB."""
        lt = np.atleast_1d(np.asarray(self.total_level, dtype=np.float64))
        return float(10.0 * np.log10(np.sum(10.0 ** (0.1 * lt))))

    def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes":
        """Plot the per-path and total normalised sound pressure levels.

        Requires matplotlib (``pip install phonometry[plot]``); returns the
        :class:`~matplotlib.axes.Axes`.
        """
        from ._plotting import plot_installed_structure_borne

        return plot_installed_structure_borne(self, ax=ax, **kwargs)


def installed_source_prediction(
    characteristic_power_level: ArrayLike,
    coupling_term: ArrayLike,
    paths: list[dict[str, Any]],
    *,
    frequencies: ArrayLike | None = None,
) -> InstalledSourceResult:
    """Predict the installed structure-borne SPL over several paths (EN 12354-5).

    :param characteristic_power_level: Characteristic level ``L_Ws,c`` (per
        band), in dB.
    :param coupling_term: Coupling term ``D_C`` (per band), in dB.
    :param paths: One dict per transmission path with keys ``adjustment_term``
        (``D_sa``), ``flanking_reduction_index`` (``R_ij,ref``) and
        ``element_area`` (``S_i``), each per band where applicable.
    :param frequencies: Band centre frequencies, in hertz, or ``None``.
    :return: The :class:`InstalledSourceResult`.
    :raises ValueError: if ``paths`` is empty.
    """
    if not paths:
        raise ValueError("'paths' must contain at least one transmission path.")
    lw_inst = installed_structure_borne_power_level(
        characteristic_power_level, coupling_term
    )
    required = ("adjustment_term", "flanking_reduction_index", "element_area")
    n_bands = np.atleast_1d(np.asarray(lw_inst, dtype=np.float64)).size
    rows = []
    for k, p in enumerate(paths):
        missing = [key for key in required if key not in p]
        if missing:
            raise ValueError(
                f"path {k} is missing required key(s): {', '.join(missing)}."
            )
        rij = np.atleast_1d(np.asarray(p["flanking_reduction_index"], dtype=np.float64))
        if rij.size not in (1, n_bands):
            raise ValueError(
                f"path {k}: 'flanking_reduction_index' has {rij.size} bands, "
                f"expected {n_bands} to match the source levels."
            )
        rows.append(
            structure_borne_pressure_level_path(
                lw_inst, p["adjustment_term"],
                p["flanking_reduction_index"], p["element_area"],
            )
        )
    path_levels = np.asarray(rows, dtype=np.float64)
    total = total_structure_borne_pressure_level(path_levels)
    freq = None if frequencies is None else np.asarray(frequencies, dtype=np.float64)
    return InstalledSourceResult(
        path_levels=path_levels,
        total_level=np.asarray(total, dtype=np.float64),
        installed_power_level=np.atleast_1d(np.asarray(lw_inst, dtype=np.float64)),
        frequencies=freq,
    )
