import ast
import re

import numpy as np
from numpy import exp, log, sqrt, sin, cos, abs as np_abs, pi
from scipy import stats
from scipy.special import ndtr  # C-level ndtr == norm.cdf w/o scipy dispatch
from scipy.stats import norm


# Restricted namespace for user-supplied repair-time expressions.  Only the
# math / distribution names actually used by the OMRAT GUI (and IWRAP-style
# func strings such as ``1 - stats.norm(loc=2, scale=1).cdf(x)``) are
# exposed; ``__builtins__`` is wiped so ``open``, ``__import__`` etc. are
# unavailable.  Combined with the AST whitelist below this stops the
# repair-function input from being a code-execution vector.
_SAFE_EVAL_GLOBALS = {
    "__builtins__": {},
    "np": np,
    "numpy": np,
    "stats": stats,
    "norm": norm,
    "exp": exp,
    "log": log,
    "sqrt": sqrt,
    "sin": sin,
    "cos": cos,
    "abs": np_abs,
    "pi": pi,
}

_SAFE_AST_NODES: tuple[type, ...] = (
    ast.Expression, ast.BinOp, ast.UnaryOp, ast.BoolOp, ast.Compare,
    ast.Constant,
    ast.Name, ast.Load,
    ast.Call, ast.Attribute, ast.keyword,
    ast.Add, ast.Sub, ast.Mult, ast.Div, ast.Mod, ast.Pow,
    ast.USub, ast.UAdd, ast.FloorDiv,
    ast.And, ast.Or, ast.Not,
    ast.Eq, ast.NotEq, ast.Lt, ast.LtE, ast.Gt, ast.GtE,
    ast.Tuple, ast.List,
)


def _validate_safe_expression(code_str: str) -> None:
    """Reject repair-time expressions that contain disallowed AST nodes."""
    tree = ast.parse(code_str, mode="eval")
    for node in ast.walk(tree):
        if not isinstance(node, _SAFE_AST_NODES):
            raise ValueError(
                f"Disallowed expression element in repair function: {type(node).__name__}"
            )


def _safe_compile(code_str: str):
    """Validate and compile a repair-time expression (raises on disallowed input)."""
    _validate_safe_expression(code_str)
    return compile(code_str, "<repair>", "eval")


def _safe_eval(code_obj, x: float) -> float:
    """Evaluate a previously :func:`_safe_compile`'d repair-time expression."""
    # Restricted namespace + AST whitelist make this eval safe.
    return eval(code_obj, _SAFE_EVAL_GLOBALS, {"x": x})  # nosec B307


def get_Fcoll(na: float, pc: float) -> float:
    """Get the accident frequncy """
    return na * pc


def get_drifting_prob(Fb: float, line_length: float, ship_speed: float) -> float:
    """Calculates the drifting probability """
    if ship_speed <= 0:
        return 1.0
    hours = line_length / ship_speed
    return 1 - exp(- Fb / (24 * 365) * hours)


def get_drift_time(distance: float, drift_speed: float) -> float:
    """Estimates the drifting time """
    if drift_speed <= 0:
        return float('inf')
    return distance / drift_speed


def repairtime_function(data, x) -> float:
    if data["active_window"] == 0:
        drift = stats.lognorm(data["std"], data["loc"], data["scale"])
        repaired = drift.cdf(x)
    else:
        repaired = _safe_eval(_safe_compile(data["func"]), x)
    return repaired


def powered_na(distance, mean_time, ship_speed):
    ai = mean_time * ship_speed
    if ai <= 0:
        return 0.0
    return exp(-distance / ai)


# Cache of (scipy dist | compiled code) keyed by a stable signature of the
# ``data`` repair-params dict.  The end-to-end profile showed ~8 s of
# ``builtins.exec`` from re-parsing ``data['func']`` on every call and
# another chunk from ``stats.lognorm(...)`` re-freezing the same
# distribution tens of thousands of times.  Caching by signature removes
# both without assuming the caller reuses the same dict object.
_REPAIR_FN_CACHE: dict[tuple, "callable"] = {}


# Recognise the two repair-function shapes that OMRAT's XML/UI writers emit.
# Matching them lets us bypass scipy's frozen-distribution dispatch (~2 ms
# per call) and use C-level ``ndtr`` / a direct Weibull formula.
_NORM_FUNC_RE = re.compile(
    r"\.norm\(loc\s*=\s*(?P<loc>-?\d+(?:\.\d+)?)\s*,\s*scale\s*=\s*(?P<scale>-?\d+(?:\.\d+)?)\s*\)\s*\.cdf\(x\)"
)
_WEIBULL_FUNC_RE = re.compile(
    r"\.weibull_min\(c\s*=\s*(?P<c>-?\d+(?:\.\d+)?)\s*,"
    r"\s*loc\s*=\s*(?P<loc>-?\d+(?:\.\d+)?)\s*,"
    r"\s*scale\s*=\s*(?P<scale>-?\d+(?:\.\d+)?)\s*\)\s*\.cdf\(x\)"
)


def _make_lognorm_repair_fn(s: float, loc: float, scale: float, drift_speed: float):
    def fn(distance: float, _s=s, _loc=loc, _scale=scale, _spd=drift_speed) -> float:
        drift_time = distance / _spd / 3600.0
        t = drift_time - _loc
        if t <= 0 or _scale <= 0:
            return 1.0
        return 1.0 - float(ndtr(log(t / _scale) / _s))
    return fn


def _make_norm_repair_fn(loc: float, scale: float, drift_speed: float):
    def fn(distance: float, _loc=loc, _scale=scale, _spd=drift_speed) -> float:
        drift_time = distance / _spd / 3600.0
        if _scale <= 0:
            return 1.0
        return 1.0 - float(ndtr((drift_time - _loc) / _scale))
    return fn


def _make_weibull_repair_fn(c: float, loc: float, scale: float, drift_speed: float):
    def fn(distance: float, _c=c, _loc=loc, _scale=scale, _spd=drift_speed) -> float:
        drift_time = distance / _spd / 3600.0
        t = drift_time - _loc
        if t <= 0 or _scale <= 0:
            return 1.0
        return float(exp(-((t / _scale) ** _c)))
    return fn


def _make_fallback_repair_fn(func_str: str, drift_speed: float):
    code = _safe_compile(func_str)

    def fn(distance: float, _code=code, _spd=drift_speed) -> float:
        x = distance / _spd / 3600.0
        return 1.0 - float(_safe_eval(_code, x))
    return fn


def _repair_fn(data: dict, drift_speed: float):
    """Return a callable ``f(distance) -> P(not repaired)`` for *data*.

    Caches the compiled scalar distribution callable, so the hot per-edge /
    per-ship call sites don't re-parse strings or re-create distribution
    objects.  For the two dominant shapes (``norm.cdf(x)`` /
    ``weibull_min.cdf(x)``) we bypass scipy's frozen-distribution machinery
    and call ``scipy.special.ndtr`` or a direct Weibull formula.
    """
    if data.get('use_lognormal') == 1:
        key: tuple = ('lognorm', float(data['std']), float(data['loc']),
                      float(data['scale']), float(drift_speed))
    else:
        key = ('func', str(data.get('func', '')), float(drift_speed))
    fn = _REPAIR_FN_CACHE.get(key)
    if fn is not None:
        return fn
    if data.get('use_lognormal') == 1:
        fn = _make_lognorm_repair_fn(float(data['std']), float(data['loc']),
                                     float(data['scale']), drift_speed)
    else:
        func_str = str(data.get('func', ''))
        m = _NORM_FUNC_RE.search(func_str)
        if m is not None:
            fn = _make_norm_repair_fn(float(m['loc']), float(m['scale']), drift_speed)
        else:
            wm = _WEIBULL_FUNC_RE.search(func_str)
            if wm is not None:
                fn = _make_weibull_repair_fn(float(wm['c']), float(wm['loc']),
                                             float(wm['scale']), drift_speed)
            else:
                fn = _make_fallback_repair_fn(func_str, drift_speed)
    _REPAIR_FN_CACHE[key] = fn
    return fn


def get_not_repaired(data: dict[str, str | float | bool], drift_speed: float, dist: float) -> float:
    """Get the probability that the ship isn't repaired."""
    return _repair_fn(data, drift_speed)(dist)


# Canonical ship-type names per OMRAT ship-type index (0-based row in the
# traffic matrix).  Used by the settings GUI and the drifting calculation to
# apply per-type parameters (e.g. blackout rate) without having to duplicate
# the mapping in multiple files.
SHIP_TYPE_NAMES: dict[int, str] = {
    0: "Other ship",
    1: "Search & rescue",
    2: "Sailing vessel",
    3: "Pleasure craft",
    4: "High speed craft",
    5: "Fishing (small)",
    6: "Fishing (large)",
    7: "Tug / push",
    8: "Passenger ferry",
    9: "Ro-pax",
    10: "Ro-ro",
    11: "Passenger cruise",
    12: "Container feeder",
    13: "Container",
    14: "General cargo",
    15: "General cargo (large)",
    16: "Dry bulk (small)",
    17: "Dry bulk (Panamax+)",
    18: "Bulk carrier",
    19: "Tanker (crude/product)",
    20: "Chemical / gas tanker",
}

# IWRAP-like default blackout rate per ship-type index, events per ship-year.
# IWRAP distinguishes "roro_passenger" (Passenger ferry, Ro-pax, Ro-ro,
# Passenger cruise) which typically have a one-order-of-magnitude lower
# blackout rate than other commercial vessels.  These match IWRAP's default
# ``blackout_other = 1.0`` and ``blackout_roro_passenger = 0.1``.
_DEFAULT_BLACKOUT_BY_SHIP_TYPE: dict[int, float] = {i: 1.0 for i in SHIP_TYPE_NAMES}
for _i in (8, 9, 10, 11):  # Passenger ferry, Ro-pax, Ro-ro, Passenger cruise
    _DEFAULT_BLACKOUT_BY_SHIP_TYPE[_i] = 0.1

# Name tokens that identify the low-blackout "roro_passenger" IWRAP class.
# Matched case-insensitively as substrings of the ship-type name so the
# defaults follow the PROJECT's actual taxonomy (e.g. the AIS type list has
# "Passenger, all ships of this type" at index 17, not 8-11).
_LOW_BLACKOUT_NAME_TOKENS: tuple[str, ...] = (
    'passenger', 'ferry', 'ro-pax', 'ropax', 'ro-ro', 'roro', 'cruise',
)


def default_blackout_by_ship_type(
    type_names: 'list[str] | None' = None,
) -> dict[int, float]:
    """Return the default per-ship-type blackout rate map (events/ship-year).

    Without arguments, returns the map for OMRAT's internal taxonomy
    (``SHIP_TYPE_NAMES``: indices 8-11 are the passenger/ro-ro classes).

    When ``type_names`` is given (the project's actual ship-type list, e.g.
    ``data['ship_categories']['types']`` -- often the AIS taxonomy where
    "Passenger" sits at index 17), the map is built by NAME matching so the
    0.1 roro_passenger rate lands on the right rows regardless of index.
    """
    if type_names is None:
        return dict(_DEFAULT_BLACKOUT_BY_SHIP_TYPE)
    out: dict[int, float] = {}
    for i, name in enumerate(type_names):
        n = str(name).lower()
        low = any(tok in n for tok in _LOW_BLACKOUT_NAME_TOKENS)
        out[i] = 0.1 if low else 1.0
    return out


# Block-coefficient defaults per OMRAT ship-type index (0-based row in traffic matrix).
# Values are typical design Cb from published statistics (IACS, literature).
# 0=Other, 1=Search & rescue, ..., 18=Cargo/bulk, 19=Tanker, 20=Passenger.
# Types with no dedicated Cb default to 0.70 (mid-range merchant vessel).
_CB_BY_SHIP_TYPE: dict[int, float] = {
    0: 0.70,   # Other ship
    1: 0.50,   # Search & rescue (fast craft)
    2: 0.65,   # Sailing vessel
    3: 0.55,   # Pleasure craft
    4: 0.65,   # High speed craft
    5: 0.60,   # Fishing (small)
    6: 0.65,   # Fishing (large)
    7: 0.65,   # Tug/push
    8: 0.70,   # Passenger ferry
    9: 0.65,   # Ro-pax
    10: 0.65,   # Ro-ro
    11: 0.70,   # Passenger cruise
    12: 0.65,   # Container feeder
    13: 0.67,   # Container
    14: 0.72,   # General cargo
    15: 0.72,   # General cargo (large)
    16: 0.72,   # Dry bulk (small)
    17: 0.80,   # Dry bulk (Panamax+)
    18: 0.80,   # Bulk carrier
    19: 0.82,   # Tanker (crude/product)
    20: 0.67,   # Chemical/gas tanker
}


def squat_m(speed_kts: float, ship_type: int = -1, cb: float | None = None) -> float:
    """Estimate under-keel squat (metres) using the Barras (1979) open-water formula.

    Barras formula (open water):
        S = C_b * V_k^2 / 100

    where:
        C_b     block coefficient (dimensionless)
        V_k     ship speed in knots

    The block coefficient is taken from ``cb`` if supplied, otherwise looked up
    by ``ship_type`` index (OMRAT row index in the traffic matrix).  If the
    ship type is unknown the default C_b = 0.75 is used (conservative).

    Returns 0 when speed <= 0.
    """
    if speed_kts <= 0.0:
        return 0.0
    if cb is None:
        cb = _CB_BY_SHIP_TYPE.get(int(ship_type) if ship_type >= 0 else -1, 0.75)
    return cb * speed_kts ** 2 / 100.0


# =============================================================================
# Ship-Ship Collision Equations (Hansen 2008 - IWRAP Theory)
# =============================================================================

def get_head_on_collision_candidates(
    Q1: float,           # Traffic in direction 1 (ships/year)
    Q2: float,           # Traffic in direction 2 (ships/year)
    V1: float,           # Speed direction 1 (m/s)
    V2: float,           # Speed direction 2 (m/s)
    mu1: float,          # Mean lateral position dir 1 (m)
    mu2: float,          # Mean lateral position dir 2 (m)
    sigma1: float,       # Std dev lateral position dir 1 (m)
    sigma2: float,       # Std dev lateral position dir 2 (m)
    B1: float,           # Beam of vessel type 1 (m)
    B2: float,           # Beam of vessel type 2 (m)
    L_w: float           # Leg segment length (m)
) -> float:
    """
    Calculate geometric number of head-on collision candidates.

    Uses Hansen Eq. 4.2-4.4:
    N_G = Q1 × Q2 × V_ij × P_G × L_w

    Where:
    - V_ij = V1 + V2 (relative closing speed for head-on)
    - P_G = Φ((μ_ij + B_ij)/σ_ij) - Φ((μ_ij - B_ij)/σ_ij)
    - μ_ij = μ1 + μ2 (mean lateral distance between vessels)
    - σ_ij = √(σ1² + σ2²) (combined standard deviation)
    - B_ij = (B1 + B2)/2 (average vessel breadth)

    Parameters
    ----------
    Q1 : float
        Traffic in direction 1 (ships/year)
    Q2 : float
        Traffic in direction 2 (ships/year)
    V1 : float
        Speed direction 1 (m/s)
    V2 : float
        Speed direction 2 (m/s)
    mu1 : float
        Mean lateral position dir 1 (m)
    mu2 : float
        Mean lateral position dir 2 (m)
    sigma1 : float
        Std dev lateral position dir 1 (m)
    sigma2 : float
        Std dev lateral position dir 2 (m)
    B1 : float
        Beam of vessel type 1 (m)
    B2 : float
        Beam of vessel type 2 (m)
    L_w : float
        Leg segment length (m)

    Returns
    -------
    float
        Geometric number of head-on collision candidates per year
    """
    # Relative closing speed for head-on collision
    V_ij = V1 + V2

    # Mean lateral separation between vessels (head-on: opposite directions).
    # Hansen Eq. 4.4: mu_ij = mu_i^(1) + mu_j^(2), both positive from own
    # sailing perspective.  OMRAT stores dir-2 means negated to a fixed
    # reference frame, so we subtract to recover the Hansen convention.
    mu_ij = mu1 - mu2

    # Combined standard deviation
    sigma_ij = sqrt(sigma1**2 + sigma2**2)

    # Average vessel breadth (collision width)
    B_ij = (B1 + B2) / 2

    # Geometric collision probability using cumulative normal distribution
    # P_G = Φ((μ_ij + B_ij)/σ_ij) - Φ((μ_ij - B_ij)/σ_ij)
    if sigma_ij > 0:
        P_G = norm.cdf((mu_ij + B_ij) / sigma_ij) - norm.cdf((mu_ij - B_ij) / sigma_ij)
    else:
        # If no variance, check if collision is certain (within beam)
        P_G = 1.0 if np_abs(mu_ij) <= B_ij else 0.0

    # Number of geometric collision candidates
    # Hansen Eq. 4.2-4.4:
    # N_G = (Q1/V1) × (Q2/V2) × V_ij × P_G × L_w
    # Where Q/V converts frequency (ships/year) to density (ships/meter)
    # This gives the correct dimension: collision candidates per year

    # Avoid division by zero
    if V1 <= 0 or V2 <= 0:
        return 0.0

    # Convert to ships per meter (density)
    # Q is ships/year, V is m/s, so Q/V gives ships/year / (m/s) = ships * s / (year * m)
    # We need to convert: ships/year / (m/year) to get ships/m
    # V_year = V * seconds_per_year
    seconds_per_year = 365.25 * 24 * 3600
    density1 = Q1 / (V1 * seconds_per_year)  # ships per meter
    density2 = Q2 / (V2 * seconds_per_year)  # ships per meter

    # N_G = density1 × density2 × V_ij × P_G × L_w × seconds_per_year
    # This gives collision candidates per year
    N_G = density1 * density2 * V_ij * P_G * L_w * seconds_per_year

    return N_G


def get_overtaking_collision_candidates(
    Q_fast: float,       # Traffic of faster vessels (ships/year)
    Q_slow: float,       # Traffic of slower vessels (ships/year)
    V_fast: float,       # Speed of faster vessel (m/s)
    V_slow: float,       # Speed of slower vessel (m/s)
    mu_fast: float,      # Mean lateral position faster (m)
    mu_slow: float,      # Mean lateral position slower (m)
    sigma_fast: float,   # Std dev lateral position faster (m)
    sigma_slow: float,   # Std dev lateral position slower (m)
    B_fast: float,       # Beam of faster vessel (m)
    B_slow: float,       # Beam of slower vessel (m)
    L_w: float           # Leg segment length (m)
) -> float:
    """
    Calculate geometric number of overtaking collision candidates.

    Similar to head-on but with:
    - V_ij = |V_fast - V_slow| (must be > 0)
    - μ_ij = μ_fast - μ_slow (same direction, subtract)

    Parameters
    ----------
    Q_fast : float
        Traffic of faster vessels (ships/year)
    Q_slow : float
        Traffic of slower vessels (ships/year)
    V_fast : float
        Speed of faster vessel (m/s)
    V_slow : float
        Speed of slower vessel (m/s)
    mu_fast : float
        Mean lateral position faster (m)
    mu_slow : float
        Mean lateral position slower (m)
    sigma_fast : float
        Std dev lateral position faster (m)
    sigma_slow : float
        Std dev lateral position slower (m)
    B_fast : float
        Beam of faster vessel (m)
    B_slow : float
        Beam of slower vessel (m)
    L_w : float
        Leg segment length (m)

    Returns
    -------
    float
        Geometric number of overtaking collision candidates per year.
        Returns 0 if V_fast <= V_slow.
    """
    # Check if overtaking is possible
    if V_fast <= V_slow:
        return 0.0

    # Relative speed for overtaking (difference in speeds)
    V_ij = V_fast - V_slow

    # Mean lateral distance (same direction: subtract)
    mu_ij = mu_fast - mu_slow

    # Combined standard deviation
    sigma_ij = sqrt(sigma_fast**2 + sigma_slow**2)

    # Average vessel breadth (collision width)
    B_ij = (B_fast + B_slow) / 2

    # Geometric collision probability
    if sigma_ij > 0:
        P_G = norm.cdf((mu_ij + B_ij) / sigma_ij) - norm.cdf((mu_ij - B_ij) / sigma_ij)
    else:
        P_G = 1.0 if np_abs(mu_ij) <= B_ij else 0.0

    # Number of geometric collision candidates
    # N_G = (Q_fast/V_fast) × (Q_slow/V_slow) × V_ij × P_G × L_w
    # Convert to density using the same approach as head-on

    # Avoid division by zero
    if V_fast <= 0 or V_slow <= 0:
        return 0.0

    seconds_per_year = 365.25 * 24 * 3600
    density_fast = Q_fast / (V_fast * seconds_per_year)  # ships per meter
    density_slow = Q_slow / (V_slow * seconds_per_year)  # ships per meter

    # N_G gives collision candidates per year
    N_G = density_fast * density_slow * V_ij * P_G * L_w * seconds_per_year

    return N_G


def get_crossing_collision_candidates(
    Q1: float,           # Traffic on leg 1 (ships/year)
    Q2: float,           # Traffic on leg 2 (ships/year)
    V1: float,           # Speed on leg 1 (m/s)
    V2: float,           # Speed on leg 2 (m/s)
    L1: float,           # Length ship type 1 (m)
    L2: float,           # Length ship type 2 (m)
    B1: float,           # Beam ship type 1 (m)
    B2: float,           # Beam ship type 2 (m)
    theta: float         # Crossing angle (radians)
) -> float:
    """
    Calculate geometric number of crossing collision candidates.

    Pedersen (1995) / Friis-Hansen (2008), crossing encounters:

    .. code-block:: text

        N_G = (Q1/V1) * (Q2/V2) * V_ij * D_ij / |sin(theta)|

    with the relative speed from the law of cosines

    .. code-block:: text

        V_ij = sqrt(V1**2 + V2**2 - 2*V1*V2*cos(theta))

    and the geometric collision diameter

    .. code-block:: text

        D_ij = (L1*V2 + L2*V1) * |sin(theta)| / V_ij
             + B1 * sqrt(1 - (V2*sin(theta)/V_ij)**2)
             + B2 * sqrt(1 - (V1*sin(theta)/V_ij)**2)

    ``D_ij`` is the width of the Minkowski sum of the two hulls
    projected onto the direction perpendicular to the *relative*
    velocity -- i.e. the swept collision cross-section seen by ship 2 in
    ship 1's frame.  Modelling each hull as an ``L x B`` rectangle, the
    projection of ship 1 onto that normal is
    ``(L1*V2*|sin th| + B1*|V1 - V2*cos th|) / V_ij``, and the identity
    ``sqrt(V_ij**2 - V2**2*sin(th)**2) == |V1 - V2*cos th|`` turns the
    beam term into the square-root form above.  Note the pairing: ``B1``
    carries ``V2`` inside the root, and vice versa.

    Sanity check -- equal speeds, ``theta = 90 deg``: every term picks up
    ``1/sqrt(2)`` and ``D_ij = (L1 + L2 + B1 + B2)/sqrt(2)``, which is
    exactly the projected width of both rectangles onto the 45-degree
    normal.

    .. note::

       Before v0.14.0 this function omitted ``V_ij`` from the numerator
       and used the unweighted ``D_ij = (L1+L2)|sin th| + (B1+B2)|cos th|``.
       That made the result scale as ``1/V**2`` instead of ``1/V`` (unlike
       head-on and overtaking, which were always correct) and dropped the
       beam contribution entirely at right-angle crossings.  Crossing and
       bend numbers from older runs are therefore not comparable with
       current ones -- they were low by roughly a factor of ``V_ij``.

    Parameters
    ----------
    Q1 : float
        Traffic on leg 1 (ships/year)
    Q2 : float
        Traffic on leg 2 (ships/year)
    V1 : float
        Speed on leg 1 (m/s)
    V2 : float
        Speed on leg 2 (m/s)
    L1 : float
        Length ship type 1 (m)
    L2 : float
        Length ship type 2 (m)
    B1 : float
        Beam ship type 1 (m)
    B2 : float
        Beam ship type 2 (m)
    theta : float
        Crossing angle (radians)

    Returns
    -------
    float
        Geometric number of crossing collision candidates per year
    """
    # Handle edge cases for crossing angle
    sin_theta = sin(theta)
    abs_sin = np_abs(sin_theta)
    if abs_sin < 1e-10:
        # Parallel or anti-parallel courses - use head-on or overtaking instead
        return 0.0

    if V1 <= 0 or V2 <= 0:
        return 0.0

    # Relative speed using law of cosines.
    V_ij = sqrt(V1**2 + V2**2 - 2 * V1 * V2 * cos(theta))
    if V_ij < 1e-10:
        return 0.0

    # Geometric collision diameter (Pedersen 1995): the width of the two
    # hulls projected onto the normal of the relative-velocity vector.
    # The lengths are weighted by the *other* ship's speed; each beam term
    # carries the other ship's speed inside the root.  ``max(0.0, ...)``
    # guards the root against round-off at theta near 0 / pi.
    d_len = (L1 * V2 + L2 * V1) * abs_sin / V_ij
    d_b1 = B1 * sqrt(max(0.0, 1.0 - (V2 * sin_theta / V_ij) ** 2))
    d_b2 = B2 * sqrt(max(0.0, 1.0 - (V1 * sin_theta / V_ij) ** 2))
    D_ij = d_len + d_b1 + d_b2

    # N_G = (Q1/V1) * (Q2/V2) * V_ij * D_ij / |sin(theta)|
    #
    # Q is ships/year and V is m/s, so ``Q / (V * seconds_per_year)`` is
    # the linear ship density in ships/m and ``V_ij * seconds_per_year``
    # is the relative speed in m/year.  The two conversions leave a single
    # ``/ seconds_per_year`` in the combined expression:
    #
    #   (ships/m) * (ships/m) * (m/year) * m  ->  ships**2/year
    seconds_per_year = 365.25 * 24 * 3600
    N_G = Q1 * Q2 * V_ij * D_ij / (V1 * V2 * abs_sin * seconds_per_year)

    return N_G


def get_bend_collision_candidates(
    Q: float,            # Traffic volume (ships/year)
    P_no_turn: float,    # Probability of not changing course at bend (≈0.01)
    L: float,            # Ship length (m)
    B: float,            # Ship beam (m)
    theta: float,        # Bend angle (radians)
    V: float,            # Ship speed (m/s)
) -> float:
    """
    Calculate bend collision candidates at waypoints.

    At bends / waypoints there is a probability that a vessel fails to
    turn and continues on the original course.  It then crosses the path
    of the traffic that did turn, at the bend's deflection angle.  This is
    the crossing encounter of two sub-streams of the *same* traffic:

    .. code-block:: text

        Q_no_turn = Q * P_no_turn          (continued straight)
        Q_turn    = Q * (1 - P_no_turn)    (made the turn)

    both at speed ``V``, crossing at ``theta``.  The result is
    :func:`get_crossing_collision_candidates` evaluated on that pair, so
    the bend inherits the Pedersen collision diameter and the ``1/V``
    scaling of the crossing formula.

    Parameters
    ----------
    Q : float
        Traffic volume (ships/year)
    P_no_turn : float
        Probability of not changing course at bend (typically ≈0.01)
    L : float
        Ship length (m)
    B : float
        Ship beam (m)
    theta : float
        Bend angle (radians) - the angle change at the waypoint
    V : float
        Ship speed (m/s).  Required: with two equal speeds the relative
        speed is ``V_ij = 2*V*|sin(theta/2)|``, which does **not** cancel
        out of the crossing formula.  Passing a placeholder ``1.0`` (as
        this function did before v0.14.0) silently drops the bend's speed
        dependence.

    Returns
    -------
    float
        Number of bend collision candidates per year
    """
    # Handle edge case of no bend
    if np_abs(theta) < 1e-10:
        return 0.0

    if V <= 0:
        return 0.0

    # Traffic that fails to turn vs traffic that turns normally.
    Q_no_turn = Q * P_no_turn
    Q_turn = Q * (1 - P_no_turn)

    # Both sub-streams are the same ship category at the same speed, so
    # the crossing formula is evaluated with L1 == L2, B1 == B2, V1 == V2.
    N_G_crossing = get_crossing_collision_candidates(
        Q1=Q_no_turn,
        Q2=Q_turn,
        V1=V,
        V2=V,
        L1=L,
        L2=L,
        B1=B,
        B2=B,
        theta=theta
    )

    return N_G_crossing


# =============================================================================
# Powered Grounding Equations (Hansen 2008, Eq. 4.15-4.17)
# =============================================================================

def get_recovery_distance(position_check_minutes: float, ship_speed: float) -> float:
    """
    Calculate the characteristic recovery distance a = λ × V.

    This is the average distance a ship travels between position checks.
    Used in powered grounding Category II calculations.

    Parameters
    ----------
    position_check_minutes : float
        Position check interval λ (minutes). Typical value is ~3 minutes.
    ship_speed : float
        Ship speed V (m/s)

    Returns
    -------
    float
        Recovery distance a (meters)
    """
    # Convert minutes to seconds for consistency with m/s speed
    position_check_seconds = position_check_minutes * 60.0
    return position_check_seconds * ship_speed


def get_powered_grounding_cat1(
    Q: float,
    Pc: float,
    prob_in_obstacle: float
) -> float:
    """
    Calculate powered grounding Category I - ships sailing directly into obstacle.

    Category I covers ships that are in the traffic lane that overlaps the obstacle.
    These ships are sailing on a course that would naturally take them into the
    obstacle zone without any course deviation.

    Based on Hansen Eq. 4.15:
    N_I = P_c × Q × prob_in_obstacle

    Where prob_in_obstacle = ∫(z_min to z_max) f_i(z) dz
    (the integral of the traffic distribution over the obstacle extent)

    Parameters
    ----------
    Q : float
        Traffic volume (ships/year)
    Pc : float
        Causation factor - probability that a ship fails to take evasive action
    prob_in_obstacle : float
        Probability of being in obstacle zone. This is the integral of the
        lateral traffic distribution f(z) over the obstacle extent from z_min
        to z_max. Typically computed as Φ((z_max - μ)/σ) - Φ((z_min - μ)/σ)
        for a normal distribution.

    Returns
    -------
    float
        Expected number of Category I powered groundings per year
    """
    return Pc * Q * prob_in_obstacle


def get_powered_grounding_cat2(
    Q: float,
    Pc: float,
    prob_at_position: float,
    distance_to_obstacle: float,
    position_check_interval: float,
    ship_speed: float
) -> float:
    """
    Calculate powered grounding Category II - ships failing to turn at bend.

    Category II covers ships that fail to change course at a bend/waypoint
    and continue on their original heading towards an obstacle. The probability
    of grounding decreases exponentially with distance from the bend as ships
    have more opportunities to detect and correct their course error.

    Based on Hansen Eq. 4.16-4.17:
    N_II = P_c × Q × f(z) × exp(-d / a)

    Where:
    - d = distance from bend to obstacle (m)
    - λ = position check interval (minutes, converted to seconds)
    - a = λ × V = distance between position checks (recovery distance)
    - f(z) = probability density at the bend position

    The exponential term exp(-d/a) represents the probability that the ship
    has NOT detected its course error before reaching the obstacle.

    Parameters
    ----------
    Q : float
        Traffic volume (ships/year)
    Pc : float
        Causation factor - probability that a ship fails to take evasive action
    prob_at_position : float
        Probability density f(z) at the bend position. This represents the
        likelihood of a ship being at the lateral position where it would
        head towards the obstacle after missing the turn.
    distance_to_obstacle : float
        Distance d from the bend to the obstacle (meters)
    position_check_interval : float
        λ - Position check interval (minutes). Typical value is ~3 minutes.
        This represents how frequently the navigator checks the ship's position.
    ship_speed : float
        Ship speed V (m/s)

    Returns
    -------
    float
        Expected number of Category II powered groundings per year
    """
    # Calculate recovery distance a = λ × V
    recovery_dist = get_recovery_distance(position_check_interval, ship_speed)

    # Avoid division by zero
    if recovery_dist <= 0:
        return 0.0

    # Exponential decay: probability of not detecting course error
    # before reaching obstacle at distance d
    prob_not_recovered = exp(-distance_to_obstacle / recovery_dist)

    return Pc * Q * prob_at_position * prob_not_recovered
