"""RBFL math helpers.

These functions implement the locked-law and detector-stage quantities used in
this reproducibility package.  They are deliberately small and transparent.
"""
from __future__ import annotations

import numpy as np


def rbfl_acceleration(g_b, A_b, a_phi):
    """Locked radial RBFL acceleration.

    Parameters
    ----------
    g_b : array-like
        Baryonic acceleration. Must be non-negative in the radial projection.
    A_b : array-like
        Baryon-derived phase operator.
    a_phi : float
        Phase-carrier acceleration scale.

    Returns
    -------
    array-like
        g_RBFL = g_b + A_b * sqrt(a_phi * g_b)
    """
    g_b = np.asarray(g_b, dtype=float)
    A_b = np.asarray(A_b, dtype=float)
    return g_b + A_b * np.sqrt(np.maximum(a_phi * g_b, 0.0))


def predicted_velocity(radius, g_rbfl):
    """Circular speed V_pred = sqrt(r g_RBFL)."""
    radius = np.asarray(radius, dtype=float)
    g_rbfl = np.asarray(g_rbfl, dtype=float)
    return np.sqrt(np.maximum(radius * g_rbfl, 0.0))


def observed_phase_amplitude(g_obs, g_b, a_phi):
    """Post-test residual phase diagnostic.

    A_phi_obs = (g_obs - g_b) / sqrt(a_phi * g_b)
    """
    g_obs = np.asarray(g_obs, dtype=float)
    g_b = np.asarray(g_b, dtype=float)
    denom = np.sqrt(np.maximum(a_phi * g_b, 1e-300))
    return (g_obs - g_b) / denom


def acceleration_ratio(g_obs, g_b):
    """Phi = g_obs / g_b."""
    g_obs = np.asarray(g_obs, dtype=float)
    g_b = np.asarray(g_b, dtype=float)
    return g_obs / np.maximum(g_b, 1e-300)


def log_detector(phi):
    """y = ln Phi."""
    phi = np.asarray(phi, dtype=float)
    return np.log(np.maximum(phi, 1e-300))


def phase_strength(sigma_b, sigma_c, p):
    """S_phase = tanh[(Sigma_b/Sigma_c)^p]."""
    sigma_b = np.asarray(sigma_b, dtype=float)
    return np.tanh(np.power(np.maximum(sigma_b / sigma_c, 0.0), p))


def phase_degree(log_sigma_gradient, z_d, r_d, r, r_phi, w_phi, eta):
    """D_phase radial projection recipe."""
    log_sigma_gradient = np.asarray(log_sigma_gradient, dtype=float)
    r = np.asarray(r, dtype=float)
    grad_term = np.exp(-eta * np.abs(log_sigma_gradient))
    geom_term = (1.0 + (z_d / r_d) ** 2) ** -0.5
    transition = (1.0 + np.exp(-(r - r_phi) / w_phi)) ** -1.0
    return grad_term * geom_term * transition


def event_band_score(lambda_phase, center=2.0, sigma=0.25):
    """Adjustable lensing/event detection-band score.

    E_event = exp[-1/2 * (ln(Lambda_phase/center)/sigma)^2]
    """
    lam = np.asarray(lambda_phase, dtype=float)
    lam = np.maximum(lam, 1e-300)
    return np.exp(-0.5 * (np.log(lam / center) / sigma) ** 2)
