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

Computation 161, the Planck spectral distribution, assembled
=============================================================
comp 160 closed four facts about light with one object, the null threshold, and
then named what was still missing: the Planck SPECTRAL distribution, which appears
nowhere in the book although its ingredients are all present.  This computation
writes that derivation down.

The assembly is ordinary statistical mechanics and is not claimed as new.  What the
computation checks is the part that is PST's: that four of the five inputs Planck's
law needs are FORCED by the substrate rather than chosen, that the fifth is inherited
rather than absent, and that the result is consistent with what the book already
prints.

THE THREE INPUTS, AND WHERE EACH COMES FROM
--------------------------------------------
  (I)   A LADDER OF COUNTABLE QUANTA, E_n = n hbar w.  INHERITED, NOT FORCED, and
        this is the one input PST does not derive.  The projected Lorentzian action
        of \S sec:projected-qft is canonically quantised in the ordinary way, once
        sublimation has supplied the timelike direction along which equal-time
        relations can be posed; the occupation number is then the eigenvalue of
        N = a^dag a, an integer.  What PST supplies is the arena, not the ladder.
        comp 160's docstring said the quantisation into countable quanta was "absent
        outright", and that was wrong: the book has it at \S sec:projected-qft, by
        the ordinary route.  That claim is corrected in comp 160 and in the body.

  (II)  BOSONIC, NOT FERMIONIC, OCCUPANCY.  The photon is an inner fluctuation of
        the gauge sector, of EVEN Clifford grade (comp 108, readings R1 and R2), so
        it is an even-order CAR composite, and even-order composites are bosonic by
        composition (\S sec:spin-statistics-structural).  The ladder therefore runs
        n = 0, 1, 2, ... and not n in {0, 1}.  Both are computed below; only the
        bosonic one gives Planck.

  (III) NO CONSERVED PHOTON NUMBER, so the chemical potential is ZERO.  This is the
        one that is easy to miss and is genuinely structural in PST.  The book keeps
        two discrete structures apart on purpose (\S sec:projected-qft): the
        TOPOLOGICAL winding number of pi_1(V) = Z, which is the conserved U(1)
        charge, and the OCCUPATION number of the Fock space, "an integer unrelated
        to the topological winding number".  Because the conserved charge is the
        winding number and NOT the occupation number, nothing in PST conserves
        photon number, so no Lagrange multiplier is attached to it and mu = 0
        exactly.  With mu != 0 the spectrum is not Planck's.

  and two carried over from comp 160:

  (IV)  LINEAR DISPERSION w = c k, from masslessness at tau_null.
  (V)   g = 2 POLARISATIONS, from the null condition k.k = 0.

WHAT IS CHECKED
---------------
  1.  mu = 0 and the bosonic ladder give <n> = 1/(e^{beta hbar w} - 1), symbolically.
      The fermionic ladder gives 1/(e^{beta hbar w} + 1) instead, and a CLASSICAL
      continuous energy gives kT, which is the ultraviolet catastrophe.
  2.  The density of states from w = ck, checked against a brute-force count of
      standing-wave modes in a cubic box.
  3.  The product is Planck's law, agreeing symbolically with the textbook forms in
      angular frequency, in frequency and in wavelength.
  4.  Integrating it over all frequencies returns EXACTLY the radiation constant
      a = g (pi^2/30) k_B^4/(hbar^3 c^3) that comp 160 checked against measurement.
      So the spectrum and the Stefan-Boltzmann coefficient are one statement at two
      resolutions, and the book's g = 2 is what makes both come out right.
  5.  Wien's displacement law in both conventions, against the CODATA constant.
  6.  The Rayleigh-Jeans and Wien limits.
  7.  The observational handle on (III): the FIRAS bound on a chemical-potential
      distortion of the CMB spectrum is |mu| < 9e-5, so PST's structural mu = 0 is
      confirmed to that precision rather than merely asserted.

A "pass": all of 1 to 7, with the integral of the assembled spectrum matching comp
160's radiation constant to machine precision and the CMB energy density at
T = 2.72548 K matching the measured value.

WHAT THIS DOES NOT DO
---------------------
The photoelectric effect is still absent.  Planck's law is a statement about the
equilibrium field; the photoelectric effect is a statement about matter absorbing one
quantum and is a different derivation, requiring the interaction vertex rather than
the free spectrum.  It remains open.
"""

import math

import sympy as sp

# CODATA 2018 / SI-exact
K_B = 1.380649e-23           # J/K            (exact)
HBAR = 1.054571817e-34       # J s            (h/2pi, h exact)
H_PLANCK = 6.62607015e-34    # J s            (exact)
C_LIGHT = 2.99792458e8       # m/s            (exact)
SIGMA_SB = 5.670374419e-8    # W m^-2 K^-4    (derived, exact)
A_RAD_MEASURED = 4.0 * SIGMA_SB / C_LIGHT
WIEN_B_CODATA = 2.897771955e-3   # m K
T_CMB = 2.72548              # K, Fixsen 2009
FIRAS_MU_BOUND = 9.0e-5      # |mu| < 9e-5, COBE/FIRAS 95% CL


def rule(char="="):
    print(char * 100)


# ----------------------------------------------------------------------------
# 1.  The occupation number, three ways
# ----------------------------------------------------------------------------
def part1_occupation():
    rule()
    print("  1.  OCCUPATION NUMBER FROM THE LADDER")
    rule()

    x = sp.Symbol("x", positive=True)      # x = beta hbar w
    q = sp.Symbol("q", positive=True)      # q = e^-x, so 0 < q < 1
    n = sp.Symbol("n", integer=True, nonnegative=True)

    # (a) bosonic ladder, mu = 0:  n = 0, 1, 2, ...  The two geometric sums are
    # taken on the q < 1 branch, which is the physical one (x > 0).
    Z_bose = sp.Sum(q ** n, (n, 0, sp.oo)).doit().args[0][0]
    M_bose = sp.Sum(n * q ** n, (n, 0, sp.oo)).doit().args[0][0]
    mean_bose = sp.simplify((M_bose / Z_bose).subs(q, sp.exp(-x)))
    target_bose = 1 / (sp.exp(x) - 1)
    ok_bose = sp.simplify(mean_bose - target_bose) == 0
    print(f"      bosonic  Z      = {Z_bose}   with q = e^-x")
    print(f"      bosonic  <n>    = {mean_bose}")
    print(f"      equals 1/(e^x - 1)                                        : {ok_bose}")

    # independent numerical check that the closed forms are the right sums
    worst_sum = 0.0
    for xv in (0.3, 1.0, 4.0):
        truncated = sum(k * math.exp(-k * xv) for k in range(0, 4000)) / sum(
            math.exp(-k * xv) for k in range(0, 4000)
        )
        closed = 1.0 / (math.exp(xv) - 1.0)
        worst_sum = max(worst_sum, abs(truncated - closed) / closed)
    ok_sum = worst_sum < 1e-12
    print(f"      truncated series agrees with the closed form (rel {worst_sum:.1e}) : {ok_sum}")

    # (b) fermionic ladder: n in {0, 1}
    Z_fermi = 1 + sp.exp(-x)
    mean_fermi = sp.simplify(sp.exp(-x) / Z_fermi)
    target_fermi = 1 / (sp.exp(x) + 1)
    ok_fermi = sp.simplify(mean_fermi - target_fermi) == 0
    print(f"      fermionic <n>   = {sp.simplify(mean_fermi)}   (would be Fermi-Dirac)")
    print(f"      equals 1/(e^x + 1)                                        : {ok_fermi}")

    # (c) classical continuous energy: equipartition, kT per mode
    E, beta = sp.symbols("E beta", positive=True)
    mean_classical = sp.simplify(
        sp.integrate(E * sp.exp(-beta * E), (E, 0, sp.oo))
        / sp.integrate(sp.exp(-beta * E), (E, 0, sp.oo))
    )
    ok_classical = sp.simplify(mean_classical - 1 / beta) == 0
    print(f"      classical <E>   = {mean_classical}   (= kT: equipartition)")
    print(f"      equals kT, which diverges when summed over modes          : {ok_classical}")

    # (d) mu != 0 changes it
    mu = sp.Symbol("mu", real=True)
    mean_mu = sp.simplify(1 / (sp.exp(x + mu) - 1))
    differs = sp.simplify(mean_mu.subs(mu, sp.Rational(1, 10)) - target_bose) != 0
    print(f"      with mu != 0    = 1/(e^(x+mu) - 1), which is NOT Planck   : {differs}")
    print()
    print("      PST forces the bosonic branch (II) and mu = 0 (III).  Neither is a choice:")
    print("      the photon is an even-order CAR composite, and the conserved U(1) charge")
    print("      is the winding number, not the occupation number, so photon number is not")
    print("      conserved and carries no Lagrange multiplier.")
    print()
    return bool(ok_bose and ok_sum and ok_fermi and ok_classical and differs)


# ----------------------------------------------------------------------------
# 2.  Density of states from the linear dispersion
# ----------------------------------------------------------------------------
def part2_density_of_states():
    rule()
    print("  2.  DENSITY OF STATES FROM w = c k")
    rule()

    # Continuum: dN/dw = g V w^2 / (2 pi^2 c^3).  Brute-force check by counting
    # the modes of a cubic box of side L with periodic boundary conditions:
    # k = (2 pi / L) (nx, ny, nz).  Count modes with |k| <= K and compare with
    # g V K^3 / (6 pi^2), the integral of the continuum density.
    g = 2
    L = 1.0
    errors = []
    for K in (40.0, 80.0, 160.0):
        n_max = int(math.ceil(K * L / (2 * math.pi))) + 1
        count = 0
        for nx in range(-n_max, n_max + 1):
            for ny in range(-n_max, n_max + 1):
                for nz in range(-n_max, n_max + 1):
                    if (2 * math.pi / L) ** 2 * (nx * nx + ny * ny + nz * nz) <= K * K:
                        count += 1
        count *= g
        predicted = g * (L ** 3) * K ** 3 / (6 * math.pi ** 2)
        rel = abs(count - predicted) / predicted
        errors.append(rel)
        print(f"      K = {K:6.1f}   counted {count:9d}   continuum {predicted:12.1f}"
              f"   rel {rel:.2e}")
    # The residual is a surface term and must shrink as the shell is refined;
    # a fixed tolerance would only be testing the largest box.
    shrinking = all(errors[i + 1] < errors[i] for i in range(len(errors) - 1))
    ok = shrinking and errors[-1] < 5e-3
    print(f"      error shrinks with K (surface term) and ends below 0.5%   : {ok}")
    print()
    return ok


# ----------------------------------------------------------------------------
# 3.  Assemble, and compare with the textbook forms
# ----------------------------------------------------------------------------
def part3_assemble():
    rule()
    print("  3.  THE ASSEMBLED SPECTRUM IS PLANCK'S LAW")
    rule()

    w, T, g = sp.symbols("omega T g", positive=True)
    hbar, kB, c, h, nu, lam = sp.symbols("hbar k_B c h nu lamda", positive=True)

    # density of states per unit volume, per unit angular frequency
    dos = g * w ** 2 / (2 * sp.pi ** 2 * c ** 3)
    # energy per mode = hbar w <n>, with <n> from part 1
    energy_per_mode = hbar * w / (sp.exp(hbar * w / (kB * T)) - 1)
    u_omega = sp.simplify(dos * energy_per_mode)
    print(f"      u(w,T) = {u_omega}")

    # textbook angular-frequency form (g = 2)
    textbook_omega = (hbar * w ** 3) / (sp.pi ** 2 * c ** 3) / (
        sp.exp(hbar * w / (kB * T)) - 1
    )
    ok_omega = sp.simplify(u_omega.subs(g, 2) - textbook_omega) == 0
    print(f"      matches the textbook u(w,T) at g = 2                      : {ok_omega}")

    # frequency form: u(nu) = u(w) dw/dnu with w = 2 pi nu, hbar = h/2pi
    u_nu = sp.simplify(
        (u_omega.subs(g, 2).subs(w, 2 * sp.pi * nu) * 2 * sp.pi).subs(
            hbar, h / (2 * sp.pi)
        )
    )
    textbook_nu = (8 * sp.pi * h * nu ** 3 / c ** 3) / (sp.exp(h * nu / (kB * T)) - 1)
    ok_nu = sp.simplify(u_nu - textbook_nu) == 0
    print(f"      u(nu,T) = {sp.simplify(u_nu)}")
    print(f"      matches 8 pi h nu^3/c^3 / (e^(h nu/kT) - 1)               : {ok_nu}")

    # wavelength form: nu = c/lam, |dnu/dlam| = c/lam^2
    u_lam = sp.simplify(textbook_nu.subs(nu, c / lam) * c / lam ** 2)
    textbook_lam = (8 * sp.pi * h * c / lam ** 5) / (sp.exp(h * c / (lam * kB * T)) - 1)
    ok_lam = sp.simplify(u_lam - textbook_lam) == 0
    print(f"      matches 8 pi h c/lam^5 / (e^(hc/lam kT) - 1)              : {ok_lam}")
    print()
    return bool(ok_omega and ok_nu and ok_lam), u_omega, (w, T, g, hbar, kB, c)


# ----------------------------------------------------------------------------
# 4.  Integrate: back to the radiation constant comp 160 checked
# ----------------------------------------------------------------------------
def part4_integrate(u_omega, syms):
    rule()
    print("  4.  INTEGRATING RETURNS comp 160's RADIATION CONSTANT")
    rule()

    w, T, g, hbar, kB, c = syms

    # sympy will not do int w^3/(e^{a w} - 1) dw directly, so use the standard
    # expansion 1/(e^y - 1) = sum_{n>=1} e^{-n y}, integrate term by term, and
    # resum.  Every step is symbolic, and the result is checked by quadrature.
    a = sp.Symbol("a", positive=True)          # a = hbar/(k_B T)
    n = sp.Symbol("n", integer=True, positive=True)
    term = sp.integrate(w ** 3 * sp.exp(-n * a * w), (w, 0, sp.oo))
    series_sum = sp.simplify(sp.summation(term, (n, 1, sp.oo)))
    print(f"      int w^3 e^(-n a w) dw = {sp.simplify(term)}")
    print(f"      summed over n >= 1    = {series_sum}   (= 6 zeta(4)/a^4)")
    ok_series = sp.simplify(series_sum - sp.pi ** 4 / (15 * a ** 4)) == 0
    print(f"      equals pi^4/(15 a^4)                                      : {ok_series}")

    # independent numerical check of that integral
    a_num = 1.0
    grid = 200000
    hi = 60.0
    step = hi / grid
    quad = sum(
        ((k + 0.5) * step) ** 3 / (math.exp(a_num * (k + 0.5) * step) - 1.0) * step
        for k in range(grid)
    )
    closed = math.pi ** 4 / (15 * a_num ** 4)
    ok_quad = abs(quad - closed) / closed < 1e-6
    print(f"      quadrature {quad:.10f}  vs  pi^4/15 = {closed:.10f}"
          f"   rel {abs(quad - closed) / closed:.2e}")
    print(f"      the series expansion is the right integral                : {ok_quad}")

    rho = sp.simplify(
        g * hbar / (2 * sp.pi ** 2 * c ** 3)
        * (sp.pi ** 4 / (15 * a ** 4)).subs(a, hbar / (kB * T))
    )
    print(f"      rho(T) = {rho}")

    target = g * sp.pi ** 2 * kB ** 4 * T ** 4 / (30 * hbar ** 3 * c ** 3)
    ok_symbolic = sp.simplify(rho - target) == 0
    print(f"      equals g (pi^2/30) k_B^4 T^4 / (hbar^3 c^3)               : {ok_symbolic}")

    # numeric, g = 2, against the measured radiation constant
    a_pred = float(
        (target / T ** 4).subs(
            {g: 2, kB: K_B, hbar: HBAR, c: C_LIGHT}
        )
    )
    rel = abs(a_pred - A_RAD_MEASURED) / A_RAD_MEASURED
    ok_numeric = rel < 1e-6
    print(f"      a(g = 2)  predicted {a_pred:.6e}  measured {A_RAD_MEASURED:.6e}"
          f"  rel {rel:.2e}")
    print(f"      reproduces comp 160's radiation constant                  : {ok_numeric}")

    a_g3 = float((target / T ** 4).subs({g: 3, kB: K_B, hbar: HBAR, c: C_LIGHT}))
    print(f"      a(g = 3) would be {a_g3:.6e}, off by a factor "
          f"{a_g3 / A_RAD_MEASURED:.3f}")
    print()
    print("      So the spectrum and the Stefan-Boltzmann coefficient are the SAME")
    print("      statement at two resolutions, and the g = 2 the null condition")
    print("      delivers is what makes both come out right.")
    print()
    return bool(ok_series and ok_quad and ok_symbolic and ok_numeric)


# ----------------------------------------------------------------------------
# 5.  Wien displacement
# ----------------------------------------------------------------------------
def part5_wien():
    rule()
    print("  5.  WIEN'S DISPLACEMENT LAW")
    rule()

    x = sp.Symbol("x", positive=True)

    # wavelength peak: d/dlam [lam^-5/(e^{a/lam} - 1)] = 0  ->  x = 5(1 - e^-x)
    x_lam = float(sp.nsolve(x - 5 * (1 - sp.exp(-x)), 5.0))
    b_pred = H_PLANCK * C_LIGHT / (x_lam * K_B)
    rel = abs(b_pred - WIEN_B_CODATA) / WIEN_B_CODATA
    ok_lam = rel < 1e-8
    print(f"      wavelength root x = {x_lam:.12f}   (5(1 - e^-x) = x)")
    print(f"      b = hc/(x k_B)    = {b_pred:.12e} m K")
    print(f"      CODATA            = {WIEN_B_CODATA:.12e} m K   rel {rel:.2e}")
    print(f"      matches CODATA                                            : {ok_lam}")

    # frequency peak: x = 3(1 - e^-x)
    x_nu = float(sp.nsolve(x - 3 * (1 - sp.exp(-x)), 3.0))
    ok_nu = abs(x_nu - 2.821439372) < 1e-8
    print(f"      frequency root  x = {x_nu:.12f}   (3(1 - e^-x) = x)")
    print(f"      the two peaks differ, as they must: they maximise different densities")
    print(f"      frequency root matches the standard value                 : {ok_nu}")
    print()
    return bool(ok_lam and ok_nu)


# ----------------------------------------------------------------------------
# 6.  The two limits
# ----------------------------------------------------------------------------
def part6_limits():
    rule()
    print("  6.  THE RAYLEIGH-JEANS AND WIEN LIMITS")
    rule()

    w, T, hbar, kB, c = sp.symbols("omega T hbar k_B c", positive=True)
    u = (hbar * w ** 3) / (sp.pi ** 2 * c ** 3) / (sp.exp(hbar * w / (kB * T)) - 1)

    # Rayleigh-Jeans: hbar w << kT.  Expand in the small parameter.
    eps = sp.Symbol("epsilon", positive=True)
    u_small = sp.simplify(sp.limit(u.subs(hbar, eps) / eps ** 0, eps, 0))
    rj = kB * T * w ** 2 / (sp.pi ** 2 * c ** 3)
    ok_rj = sp.simplify(u_small - rj) == 0
    print(f"      hbar w << kT  ->  {sp.simplify(u_small)}")
    print(f"      equals k_B T w^2/(pi^2 c^3), the Rayleigh-Jeans law        : {ok_rj}")
    print("      and that is the ultraviolet catastrophe: w^2 is not integrable to infinity.")

    # Wien tail: hbar w >> kT
    x = sp.Symbol("x", positive=True)
    ratio = sp.limit(
        (1 / (sp.exp(x) - 1)) / sp.exp(-x), x, sp.oo
    )
    ok_wien = sp.simplify(ratio - 1) == 0
    print(f"      hbar w >> kT  ->  occupation falls as e^-x, ratio -> {ratio}")
    print(f"      recovers Wien's exponential tail                          : {ok_wien}")
    print()
    return bool(ok_rj and ok_wien)


# ----------------------------------------------------------------------------
# 7.  The CMB, and the observational handle on mu = 0
# ----------------------------------------------------------------------------
def part7_cmb():
    rule()
    print("  7.  THE CMB: mu = 0 IS MEASURED, NOT ONLY ASSERTED")
    rule()

    rho_pred = A_RAD_MEASURED * T_CMB ** 4
    print(f"      T_CMB = {T_CMB} K  ->  rho = a T^4 = {rho_pred:.6e} J m^-3")
    print(f"      measured CMB energy density is 4.175e-14 J m^-3")
    ok_rho = abs(rho_pred - 4.175e-14) / 4.175e-14 < 2e-3
    print(f"      agrees with the measured value                            : {ok_rho}")

    # What a nonzero mu would do, at the FIRAS bound, in the Wien tail where the
    # distortion is largest relative to the spectrum.
    def occ(x, mu):
        return 1.0 / (math.exp(x + mu) - 1.0)

    print()
    print("      a chemical-potential distortion n = 1/(e^(x+mu) - 1), at the FIRAS bound:")
    worst = 0.0
    for x in (0.5, 1.0, 3.0, 10.0):
        frac = abs(occ(x, FIRAS_MU_BOUND) - occ(x, 0.0)) / occ(x, 0.0)
        worst = max(worst, frac)
        print(f"        x = {x:5.1f}   fractional change at |mu| = {FIRAS_MU_BOUND:.0e}"
              f"   : {frac:.3e}")
    ok_mu = worst < 1e-3
    print(f"      COBE/FIRAS bounds |mu| < {FIRAS_MU_BOUND:.0e}, so the spectrum is Planck")
    print(f"      to better than a part in 10^4.  PST does not fit mu; it has no")
    print(f"      conserved photon number to attach one to, so mu = 0 structurally.")
    print(f"      structural prediction confirmed within the bound           : {ok_mu}")
    print()
    return bool(ok_rho and ok_mu)


def main():
    print()
    rule()
    print("  COMPUTATION 161: THE PLANCK SPECTRAL DISTRIBUTION, ASSEMBLED")
    rule()
    print("  The assembly is ordinary statistical mechanics.  What is PST's is that four")
    print("  of the five inputs are forced: the bosonic branch, mu = 0, the linear")
    print("  dispersion and g = 2.  The countable ladder is inherited, not derived.")
    print()

    c1 = part1_occupation()
    c2 = part2_density_of_states()
    c3, u_omega, syms = part3_assemble()
    c4 = part4_integrate(u_omega, syms)
    c5 = part5_wien()
    c6 = part6_limits()
    c7 = part7_cmb()

    rule()
    print("  ASSESSMENT")
    rule()
    print(f"  1. occupation number: bosonic, fermionic and classical branches  : {c1}")
    print(f"  2. density of states from w = ck, against a lattice count        : {c2}")
    print(f"  3. the product is Planck's law in all three conventions          : {c3}")
    print(f"  4. integrating returns comp 160's radiation constant exactly     : {c4}")
    print(f"  5. Wien displacement, both conventions, against CODATA           : {c5}")
    print(f"  6. Rayleigh-Jeans and Wien limits                                : {c6}")
    print(f"  7. the CMB, and the FIRAS bound on mu                            : {c7}")
    print()
    ok = all([c1, c2, c3, c4, c5, c6, c7])
    print("  RESULT: " + (
        "THE PLANCK SPECTRUM ASSEMBLES FROM FIVE INPUTS, FOUR OF THEM FORCED BY PST -- "
        "the bosonic branch from even CAR-composition order, mu = 0 from the winding "
        "number rather than the occupation number being the conserved charge, the linear "
        "dispersion from masslessness at tau_null, and g = 2 from the null condition; "
        "the countable ladder is inherited from canonical quantisation of the projected "
        "action rather than derived.  Integrating the result returns the radiation "
        "constant the book already prints"
        if ok else "MISMATCH -- revise"))
    print()


if __name__ == "__main__":
    main()