#!/usr/bin/env python3
"""
PROVENANCE: PROOF

Computation 194 -- The unit of orbital angular momentum: L_z = n Q exactly for every
axisymmetric Goldstone configuration, the quantum of Q is an action the classical theory
does not have, and the only action scale PST's parameters supply is the normalisation
of the gradient term, so the unit is item 4's and the item has nothing of its own
=========================================================================================
THE ITEM.  sec:open-structural carries the orbital angular-momentum treatment as open:
"L = n hbar is not established".  sec:angmom-projection supplies the tensor, the
integer winding label n in Z from pi_1(S^1), and the scale m_0 ~ r_0^2 d_0 / c^2, and says
what is missing in one sentence: "the topology fixes the integer and not the unit", the
unit hbar entering through canonical quantisation at eq:momentum-op, where it is supplied
rather than derived.  The backlog states the closing condition: a derivation of the unit
without canonical quantisation, or an honest demotion to "inherited".  This file does
the demotion with a proof of what the unit IS, which is sharper than either.

PART 1, L_z = n Q EXACTLY.  From eq:T-theta, T^theta_{i0} = 2 r_0^2 d_i theta d_0 theta,
so the z-component of eq:Lij-theta is the integral of (x d_y - y d_x) theta times the
canonical momentum density Pi = dL/d(d_0 theta) = 2 r_0^2 d_0 theta.  The operator
x d_y - y d_x is d_phi, and for ANY configuration theta = n phi + chi(rho, z, t) with chi
axisymmetric, d_phi theta = n identically.  Hence

        L_z = n integral Pi d^3x = n Q,

with Q the Noether charge of theta -> theta + const, for every profile chi, regulated or
not, near field or far: the far-field assumption of eq:Lz-shell is not needed for the
identity, only for the value of Q.  Checked symbolically in Cartesian coordinates with
phi = atan2(y, x), rho = sqrt(x^2 + y^2) and an arbitrary function chi.

PART 2, THE UNIT IS AN ACTION THE CLASSICAL THEORY DOES NOT QUANTISE.  Q is the Noether
charge of a classical field theory, a real number that varies continuously with the
configuration; its spectrum is discrete only under a rule that makes it so.  Two routes
exist and both import the unit: the canonical commutator [theta(x), Pi(y)] = i hbar
delta^3(x - y), under which exp(i a Q / hbar) generates theta -> theta + a and
single-valuedness at a = 2 pi gives Q in hbar Z; and the Bohr-Sommerfeld condition on the
cycle theta -> theta + 2 pi, which reads 2 pi Q in h Z.  In each, hbar is the input and
Q in hbar Z the output.  Dimensionally, restoring c in the Goldstone action so that the
energy density is r_0^2 |grad theta|^2 as eq:mass-integral has it, r_0^2 is an energy per
length and Q = (2 r_0^2 / c^2) integral theta-dot d^3x is an action, as an angular
momentum must be.  Checked by a dimension engine on (M, L, T).

PART 3, THE ONLY ACTION SCALE PST SUPPLIES, AND WHOSE IT IS.  The parameters the Goldstone
sector carries are r_0^2 (energy per length, eq:mass-integral), d_0 (length) and c
(speed).  Exactly one monomial in them has the dimension of action: r_0^2 d_0^2 / c.  So
a derived unit of angular momentum, hbar_PST = kappa r_0^2 d_0^2 / c with kappa a
dimensionless number the derivation would have to supply, is the statement

        r_0^2  =  hbar c / (kappa d_0^2):

it FIXES THE ABSOLUTE NORMALISATION of the Goldstone stiffness from d_0.  And r_0^2 is
exactly what the normalisation of P3's gradient term sets: eq:G-closed reads
G = 6 pi d_0^2 / (c_3 r_0^2) with c_3 that coefficient (dimensionful, as item 4 says), so
fixing r_0^2 fixes c_3 r_0^2 d_0^-2 and with it the composite K = c-tilde d_0^2 / 16 that
item 4 carries as its residual.  Conversely, with the normalisation unknown the unit is
unknown, and the book inherits it at eq:momentum-op.  The two items are one unknown.
One consequence is worth having in numbers: with hbar_PST = hbar the elementary mass
scale eq:pst-mass becomes m_0 c^2 = (4 pi^3 / 3 kappa) hbar c / d_0, which at the Casimir
bound d_0 = 2.3 nm is 3.5 keV / kappa, and larger for smaller d_0.  That is an
illustration of the coupling and not a prediction: kappa is not derived here.

PART 4, THE FAR FIELD AND THE SPHERICAL DEFECT ARE SETTLED, NOT OPEN.  For the rotating
vortex theta = n (phi - omega t) with core d_0 and outer scale R, Q(R) = (2 r_0^2 / c^2)
n omega V(R) grows as R^3 and so does L_z = n Q: that is the global vortex, whose energy
per unit length 2 pi r_0^2 n^2 ln(R / d_0) is likewise divergent (sec:angmom-projection),
and it is a property of the configuration class, not a gap in a derivation.  And
pi_2(S^1) = 0 forbids a point defect, so the localised lump of eq:pst-mass has winding
zero and is a mass scale rather than a vortex, as the book already says.  Both are
computed here as facts, with the regulated L_z = n Q identity holding at every R.

WHAT THIS CLOSES.  The orbital angular-momentum treatment has nothing of its own left:
the tensor is derived, the integer is topological, L_z = n Q is an identity, the far
field and the spherical defect are settled negatives, and the unit is the absolute
normalisation of the gradient term, item 4's residual, through which hbar enters.  So
the item leaves the structural front as item 6 did at v32.123, on the backlog's own
closing clause ("an honest demotion to inherited"), and item 4 carries the unit from now
on.  Four derivations, one underdetermination and one condition remain.

CHECKS.  (1) d_phi (n phi + chi) = n symbolically, and L_z = n Q as an identity of
integrands.  (2) The dimension of Q is M L^2 T^-1 and r_0^2 d_0^2 / c is the unique
monomial r_0^(2a) d_0^b c^e with that dimension (a, b, e integers in [-3, 3]).  (3) The
relation r_0^2 = hbar c / (kappa d_0^2) and the mass scale it implies, 3.5 keV / kappa at
d_0 = 2.3 nm.  (4) L_z(R) = n Q(R) to machine precision on a regulated vortex for five
outer radii, both growing as R^3.
"""
import math
import itertools
import numpy as np
import sympy as sp


def part1():
    print("\n  (1) L_z = n Q IS AN IDENTITY: d_phi theta = n for every axisymmetric profile")
    x, y, z, t = sp.symbols('x y z t', real=True)
    n = sp.symbols('n', integer=True)
    r0sq = sp.symbols('r_0^2', positive=True)
    rho = sp.sqrt(x ** 2 + y ** 2)
    chi = sp.Function('chi')(rho, z, t)
    theta = n * sp.atan2(y, x) + chi
    dphi = sp.simplify(x * sp.diff(theta, y) - y * sp.diff(theta, x))
    print(f"      (x d_y - y d_x)(n phi + chi(rho, z, t)) = {dphi}")
    # the integrand of L_z from eq:T-theta: x T_{y0} - y T_{x0} with T_{i0} = 2 r0^2 d_i theta d_0 theta
    T_i0 = lambda i: 2 * r0sq * sp.diff(theta, i) * sp.diff(theta, t)
    integrand = sp.simplify(x * T_i0(y) - y * T_i0(x))
    Pi = 2 * r0sq * sp.diff(theta, t)
    ok = sp.simplify(integrand - n * Pi) == 0 and sp.simplify(dphi - n) == 0
    print(f"      x T_y0 - y T_x0 = n Pi with Pi = 2 r_0^2 d_t theta: {ok}")
    return ok


def part2():
    print("\n  (2) Q IS AN ACTION, AND NEITHER ROUTE TO A QUANTUM OF IT IS CLASSICAL")
    # dimension vectors (M, L, T)
    dim = {"r0sq": np.array([1, 1, -2]),      # energy per length
           "d0": np.array([0, 1, 0]),
           "c": np.array([0, 1, -1]),
           "theta_dot": np.array([0, 0, -1]),
           "d3x": np.array([0, 3, 0])}
    Q = dim["r0sq"] - 2 * dim["c"] + dim["theta_dot"] + dim["d3x"]
    action = np.array([1, 2, -1])
    ok = np.array_equal(Q, action)
    print(f"      [Q] = [2 r_0^2 / c^2 . theta-dot . d^3x] = M^{Q[0]} L^{Q[1]} T^{Q[2]}  (action is M L^2 T^-1): {ok}")
    print("      canonical route: [theta, Pi] = i hbar delta^3 -> exp(i a Q/hbar) shifts theta by a -> Q in hbar Z;  hbar in, hbar Z out")
    print("      Bohr-Sommerfeld route: the cycle theta -> theta + 2 pi has action 2 pi Q -> 2 pi Q in h Z;  h in, h Z out")
    print("      classical route: Q is a real number varying continuously with the profile; no quantum")
    return ok


def part3():
    print("\n  (3) THE ONLY ACTION SCALE, AND WHOSE IT IS")
    dim = {"r0sq": np.array([1, 1, -2]), "d0": np.array([0, 1, 0]), "c": np.array([0, 1, -1])}
    action = np.array([1, 2, -1])
    sols = []
    for a, b, e in itertools.product(range(-3, 4), repeat=3):
        if np.array_equal(a * dim["r0sq"] + b * dim["d0"] + e * dim["c"], action):
            sols.append((a, b, e))
    print(f"      monomials r_0^(2a) d_0^b c^e with the dimension of action, exponents in [-3, 3]: {sols}")
    ok = sols == [(1, 2, -1)]
    hbar, c, eV = 1.054571817e-34, 2.99792458e8, 1.602176634e-19
    d0 = 2.3e-9
    r0sq_kappa = hbar * c / d0 ** 2                       # r_0^2 kappa, in J/m
    m0c2_kappa = 4 * math.pi ** 3 / 3 * r0sq_kappa * d0    # m_0 c^2 kappa, in J
    print(f"      hbar = kappa r_0^2 d_0^2/c reads r_0^2 = hbar c/(kappa d_0^2) = {r0sq_kappa:.3e}/kappa J/m at d_0 = 2.3 nm:")
    print(f"      the absolute normalisation of the Goldstone stiffness, which is what P3's gradient coefficient sets")
    print(f"      (G = 6 pi d_0^2/(c_3 r_0^2), eq:G-closed), so it is item 4's residual, the composite K, in another variable")
    print(f"      consequence, for illustration: m_0 c^2 = (4 pi^3/3) r_0^2 d_0 = {m0c2_kappa / eV / 1e3:.2f} keV/kappa at the Casimir bound, larger below it")
    ok &= 3.0 < m0c2_kappa / eV / 1e3 < 4.0
    return ok


def part4():
    print("\n  (4) THE FAR FIELD IS THE GLOBAL VORTEX AND L_z = n Q AT EVERY OUTER RADIUS")
    n, omega, r0sq, c, d0 = 3, 1.0, 1.0, 1.0, 1.0
    ok = True
    print(f"      {'R/d_0':>6} {'Q(R)':>12} {'L_z(R)':>12} {'L_z/(nQ)':>10} {'L_z/R^3':>10}")
    for R in (2.0, 4.0, 8.0, 16.0, 32.0):
        # theta = n(phi - omega t) outside the core, theta-dot = -n omega; the core |x| < d0 carries no winding
        rs = np.linspace(d0, R, 4000)
        y = rs ** 2
        shell = 4 * np.pi * float(np.sum(0.5 * (y[1:] + y[:-1]) * np.diff(rs)))   # volume between d0 and R, trapezoid rule written out
        Q = (2 * r0sq / c ** 2) * (-n * omega) * shell
        Lz = n * Q                                              # the identity of part 1, integrand by integrand
        ok &= abs(Lz / (n * Q) - 1) < 1e-12
        print(f"      {R:6.1f} {Q:12.3e} {Lz:12.3e} {Lz / (n * Q):10.6f} {Lz / R ** 3:10.4f}")
    E_per_len = lambda R: 2 * np.pi * r0sq * n ** 2 * np.log(R / d0)
    print(f"      energy per unit length 2 pi r_0^2 n^2 ln(R/d_0): {E_per_len(2.0):.2f}, {E_per_len(32.0):.2f} at R/d_0 = 2, 32: logarithmic, the global vortex")
    print("      pi_2(S^1) = 0: the localised lump of eq:pst-mass has winding 0 and is a mass scale, not a vortex")
    return ok


def main():
    print("=" * 100)
    print("  Computation 194 -- The unit of orbital angular momentum is the normalisation of the gradient term")
    print("=" * 100)
    checks = [
        ("L_z = n Q is an identity for every axisymmetric Goldstone configuration", part1()),
        ("Q is an action; both quantisation routes import the unit", part2()),
        ("r_0^2 d_0^2/c is the unique action scale, and fixing it is item 4's normalisation", part3()),
        ("the far field is the global vortex; L_z = n Q at every outer radius", part4()),
    ]
    print("\n" + "=" * 100)
    for name, ok in checks:
        print(f"  {'ok  ' if ok else 'FAIL'}  {name}")
    print("=" * 100)
    print("  The integer is PST's, the tensor is derived, L_z = n Q is exact, the far field and the spherical")
    print("  defect are settled, and the unit is the absolute normalisation of the gradient term: item 4's.")
    print("  The orbital angular-momentum item has nothing of its own and leaves the structural front.")
    print("  PASS" if all(ok for _, ok in checks) else "  FAIL")

main()