import mpmath as mp
mp.mp.dps=30
def n(r): return mp.e**(2/r)
def f(r): return n(r)*r
def r0(b): return mp.findroot(lambda r: f(r)-b, max(b,mp.mpf('2.0001')))
def alpha(b):
    rr0=r0(b)
    # substitute r = rr0/(1-z), z in [0,1), dr = rr0/(1-z)^2 dz, r->inf as z->1
    def g(z):
        r=rr0/(1-z)
        dr=rr0/(1-z)**2
        denom=f(r)**2-b**2
        if denom<=0: denom=mp.mpf('1e-50')
        return b/(r*mp.sqrt(denom))*dr
    return 2*mp.quad(g,[0,1])-mp.pi
bc=2*mp.e
for db in ['1e-3','1e-5','1e-7','1e-9']:
    x=mp.mpf(db); b=bc*(1+x)
    print(f"x={db}: alpha={mp.nstr(alpha(b),14)}  alpha+ln(x)={mp.nstr(alpha(b)+mp.log(x),14)}")
# fit
x1,x2=mp.mpf('1e-9'),mp.mpf('1e-11')
a1=alpha(bc*(1+x1)); a2=alpha(bc*(1+x2))
abar=(a2-a1)/(mp.log(x1)-mp.log(x2))
bbar=a1+abar*mp.log(x1)
print("\nabar=",mp.nstr(abar,14)," bbar=",mp.nstr(bbar,14))
print("gamma_L (Lyapunov per orbit, =2pi/?)...")
# analytic Lyapunov per phi for unstable photon orbit:
rph=mp.mpf(2)
lam = mp.sqrt(mp.diff(f,rph,2)*rph**2/(2*f(rph)))
print("lam_phi=",mp.nstr(lam,14)," 1/lam_phi=",mp.nstr(1/lam,14))
