"""Own razor-thin exponential benchmark; never substituted for native SPARC disks."""
import math
from scipy.integrate import quad
from scipy.special import ive,kve,ellipkm1


def speed_squared_dimensionless(radius):
    if not math.isfinite(radius) or radius<=0:raise ValueError('Positive dimensionless radius required')
    y=radius/2
    return 2*y*y*(ive(0,y)*kve(0,y)-ive(1,y)*kve(1,y))

def ring_potential(radius):
    """Phi/(GM/Rd), directly integrating Newtonian circular-ring potential."""
    if radius<=0:raise ValueError('Positive radius required')
    def integrand(source_radius):
        complement=((source_radius-radius)/(source_radius+radius))**2
        return math.exp(-source_radius)*source_radius/(source_radius+radius)*ellipkm1(complement)
    # The logarithmic ring singularity is an integrable endpoint on both pieces.
    inner=quad(integrand,0,radius,epsabs=2e-13,epsrel=2e-13,limit=300)[0]
    outer=quad(integrand,radius,60,epsabs=2e-13,epsrel=2e-13,limit=300)[0]
    return -2/math.pi*(inner+outer)

def ring_force_speed_squared(radius,relative_step):
    step=relative_step*max(radius,.1)
    derivative=(-ring_potential(radius+2*step)+8*ring_potential(radius+step)-8*ring_potential(radius-step)+ring_potential(radius-2*step))/(12*step)
    return radius*derivative