"""Independent dimensional conventions and 3-D Roche-lobe volume quadrature.

G = total mass = separation = orbital angular frequency = 1. The donor is
at x=-mu and q=donor/accretor. Positive Omega is minus the Roche potential.
No centrifugal z**2 term is present. Rays stop at the first donor-connected
critical surface; the L1 cusp itself has zero solid-angle measure.
"""
from __future__ import annotations

import numpy as np
from functools import lru_cache
from numpy.polynomial.legendre import leggauss
from scipy.optimize import brentq


def omega(mu, x, y, z=0):
    return ((x*x+y*y)/2 + (1-mu)/np.sqrt((x+mu)**2+y*y+z*z)
            + mu/np.sqrt((x-1+mu)**2+y*y+z*z))


def l1(q: float) -> float:
    mu = 1/(1+q)
    def acceleration(x):
        return x-(1-mu)*(x+mu)/abs(x+mu)**3-mu*(x-1+mu)/abs(x-1+mu)**3
    return brentq(acceleration, -mu+1e-12, 1-mu-1e-12, xtol=1e-14)


def eggleton(q):
    cube_root = np.cbrt(q)
    return .49*cube_root*cube_root/(.6*cube_root*cube_root+np.log1p(cube_root))


def radius_log_derivative(q):
    cube_root = np.cbrt(q)
    return 2/3-(1.2*cube_root*cube_root+cube_root/(1+cube_root))/(3*(.6*cube_root*cube_root+np.log1p(cube_root)))


@lru_cache(maxsize=512)
def volume_response(q: float, order: int) -> tuple[float, float]:
    """Radius and dln(radius)/dln(q), differentiating the implicit surface.

    The critical potential's derivative has no dL1 term because its radial
    gradient vanishes there. This response is not the Eggleton fit derivative.
    Nodes exclude the zero-measure L1 cusp; angular refinement is still required.
    """
    mu = 1/(1+q)
    critical_x = l1(q)
    critical_omega = omega(mu, critical_x, 0)
    polar_cosine, polar_weights = leggauss(order)
    azimuth_nodes, azimuth_weights = leggauss(order)
    azimuth = (azimuth_nodes+1)*np.pi/4
    transverse = np.sqrt(1-polar_cosine[:, None]**2)
    directions = (polar_cosine[:, None], transverse*np.cos(azimuth), transverse*np.sin(azimuth))
    upper = np.full((order, order), critical_x+mu)
    lower = np.zeros_like(upper)
    def on_ray(radius):
        return omega(mu, -mu+radius*directions[0], radius*directions[1], radius*directions[2])
    if np.any(on_ray(upper) > critical_omega+1e-12):
        raise ValueError("L1 distance fails to bracket the donor-connected surface")
    for _ in range(58):
        middle = (lower+upper)/2
        interior = on_ray(middle) > critical_omega
        lower = np.where(interior, middle, lower)
        upper = np.where(interior, upper, middle)
    radius = (upper+lower)/2
    weights = polar_weights[:, None]*azimuth_weights[None, :]
    radial_cube_sum = np.sum(radius**3*weights)
    donor_l1_distance = critical_x+mu
    companion_distance = np.sqrt(1+radius**2-2*radius*directions[0])
    potential_mu = -1/radius+1/companion_distance+mu-radius*directions[0]
    critical_mu = -1/donor_l1_distance+1/(1-donor_l1_distance)+mu-donor_l1_distance
    potential_radial = (-(1-mu)/radius**2-mu*(radius-directions[0])/companion_distance**3
                        -mu*directions[0]+radius*(directions[0]**2+directions[1]**2))
    radial_logq_derivative = -mu*(1-mu)*(critical_mu-potential_mu)/potential_radial
    logarithmic_response = np.sum(radius**2*radial_logq_derivative*weights)/radial_cube_sum
    return float((radial_cube_sum/4)**(1/3)), float(logarithmic_response)


def volume_radius(q: float, order: int) -> float:
    return volume_response(q, order)[0]