"""A rotated whole-sphere quadrature used only as an independent reference.

Following the separate AC22 audit, choose the z polar axis and midpoint periodic
azimuth, without the primary calculation's x-axis/octant factorization. Finite
differences of these volumes independently test the primary implicit derivative.
No code or result is imported from /tmp or another worktree.
"""
from functools import lru_cache
import numpy as np
from numpy.polynomial.legendre import leggauss
from scipy.optimize import brentq


@lru_cache(maxsize=512)
def independent_radius(q, order):
    donor_fraction = q/(1+q)
    companion_fraction = 1/(1+q)
    def potential(x, y, z):
        return (donor_fraction/np.sqrt((x+companion_fraction)**2+y*y+z*z)
                + companion_fraction/np.sqrt((x-donor_fraction)**2+y*y+z*z)
                + (x*x+y*y)/2)
    def acceleration(x):
        donor_offset = x+companion_fraction
        companion_offset = x-donor_fraction
        return x-donor_fraction*donor_offset/abs(donor_offset)**3-companion_fraction*companion_offset/abs(companion_offset)**3
    critical_x = brentq(acceleration, -companion_fraction+1e-11, donor_fraction-1e-11, xtol=1e-14)
    critical = potential(critical_x, 0, 0)
    z_cosine, weights = leggauss(order)
    azimuth = (np.arange(2*order)+.5)*np.pi/order
    transverse = np.sqrt(1-z_cosine*z_cosine)[:, None]
    direction_x, direction_y = transverse*np.cos(azimuth), transverse*np.sin(azimuth)
    direction_z = z_cosine[:, None]
    lower = np.zeros_like(direction_x)
    upper = np.full_like(direction_x, critical_x+companion_fraction)
    if np.max(potential(-companion_fraction+upper*direction_x, upper*direction_y, upper*direction_z)-critical) > 1e-11:
        raise ValueError("Independent donor-connected surface was not bracketed")
    for _ in range(54):
        middle = (lower+upper)/2
        interior = potential(-companion_fraction+middle*direction_x, middle*direction_y, middle*direction_z) > critical
        lower = np.where(interior, middle, lower)
        upper = np.where(interior, upper, middle)
    volume_radius_cubed = np.sum(weights[:, None]*((lower+upper)/2)**3)/(4*order)
    return float(volume_radius_cubed**(1/3))


def independent_log_derivative(q, order, log_step):
    return (np.log(independent_radius(q*np.exp(log_step), order))
            - np.log(independent_radius(q*np.exp(-log_step), order)))/(2*log_step)