"""Check of the polar density in l:sigma-radial-integral (Section 6).

Claim.  In the polar coordinates (e:sigma-polar) of the rank-one target space,
every generator field Yhat, Y = E_ij in gl(2|2), built from (e:sigma-action) by
casimir.py, is divergence free for the density rho = d^{-2}, d = B - F:

    div_rho(D) = - sum_{x in (B, F, theta_B, theta_F)} d_x(rho v^x)
                 + sum_a d_{xi_a}(rho v^{xi_a})  = 0,

where D = sum_I v^I d_I is the even field theta*Yhat (odd Y) or Yhat (even Y),
and d_{xi_a} is the left derivative.  This is the condition for
int d^{-2} [D f]_top dB dF dtheta_B dtheta_F = 0 for all compactly supported f.

Part 1 checks the identity for rho = d^{-2} in exact arithmetic (all 16
generators, every Grassmann coefficient).  Part 2 evaluates the divergence for
rho = d^{-m}, m = 1, 2, 3, at two sample points, to show that the exponent 2 is
the only one that works.

Run:  python3 check_polar_density.py      (from any directory)
"""
import sys, os
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import sympy as sp
from grass import G
from casimir import basis, B, F, tB, tF, odd_gens, even_coords


def divergence(v, rho):
    div = G()
    for I in range(4):
        div = div - (v[I] * rho).diff(even_coords[I])
    for a in range(4):
        div = div + (v[4 + a] * rho).dleft(odd_gens[a])
    return div


# Part 1: exact check for d^{-2}
rho = (B - F) ** (-2)
bad = 0
for (ij, p, v) in basis:
    for mono, c in divergence(v, rho).d.items():
        s = sp.simplify(sp.expand(c))
        if s != 0 and sp.simplify(s.rewrite(sp.exp)) != 0:
            bad += 1
            print("nonzero:", ij, mono, s)
print("Part 1 (exact, d^-2): generators:", len(basis),
      " nonzero divergence coefficients:", bad)

# Part 2: numerical comparison of the exponents 1, 2, 3
pts = [{B: sp.Rational(23, 10), F: sp.Rational(3, 10), tB: sp.Rational(7, 10), tF: sp.Rational(-4, 10)},
       {B: sp.Rational(51, 10), F: sp.Rational(-7, 10), tB: sp.Rational(1, 3), tF: sp.Rational(2, 1)}]
for m in [1, 2, 3]:
    rho = (B - F) ** (-m)
    worst = 0
    for (ij, p, v) in basis:
        div = divergence(v, rho)
        for pt in pts:
            for mono, c in div.d.items():
                worst = max(worst, abs(complex(sp.N(c.subs(pt), 30))))
    print(f"Part 2: density d^-{m}: max |div| over the 16 generators and 2 points = {worst:.3e}")