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

Computation 160, the photon: what tau_null fixes, and what is still missing
===========================================================================
The photon appears in the book as the U(1) connection (\S sec:gauge), as exactly
massless, and as the thing that runs along tau_null (\S sec:black-holes).  The
particle side is thinner: polarisation is mentioned three times, helicity once,
blackbody once, the photoelectric effect not at all.  comp 159 then found that a
photon contributes nothing to the stress tensor from its own worldline, because the
derived weight vanishes on null separations.  This computation follows that thread and
finds that ONE object, the null threshold, fixes four separate things about light.

WHAT tau_null FIXES
-------------------
  (1) THE LIGHT CONE AND THE SIGNATURE.  eq:metric: a pair with
      T({a,b}) = tau_null is null.  This is the book's, not new here.

  (2) EXACTLY TWO POLARISATIONS.  The physical states of a vector excitation with
      momentum k are {eps : k.eps = 0} modulo eps ~ eps + lambda k.  The transverse
      subspace is three-dimensional for any k.  The gauge direction can be quotiented
      out ONLY IF k lies in its own transverse subspace, which is k.k = 0 exactly.
      So a null k gives 3 - 1 = 2 states and a timelike k gives 3 - 0 = 3, verified
      over 200 random momenta of each kind.  The discriminant is the null condition,
      which is the one tau_null defines.  The counting itself is ordinary gauge
      theory; what PST supplies is the input, that the photon's momentum is null,
      via the unbroken U(1) of \S sec:gauge.

  (3) THE EQUATION OF STATE, w = 1/3.  comp 159: a null worldline has s = 0 along it,
      so the derived weight W(s) = -|s| of comp 158 vanishes identically on it.  The
      radiation value is forced by the weight, not obtained as a limit.

  (4) THE STEFAN-BOLTZMANN LAW.  With p = rho/3 the first law is integrable only if
      T is proportional to rho^(1/4), so rho is proportional to T^4.  This is
      Boltzmann's 1884 argument and it needs no quantum input at all; solved
      symbolically below.  So PST reaches Stefan-Boltzmann from (3), which it derives,
      rather than assuming it.

  (5) AND THE TWO OF (2) FIXES THE CONSTANT.  The radiation constant is
      a = g (pi^2/30) k_B^4 / (hbar^3 c^3) with g the number of polarisations.  With
      g = 2 it reproduces the measured 7.5657e-16 J m^-3 K^-4; with g = 3 it would be
      half as large again.  So the null condition fixes not only the T^4 scaling but
      its coefficient.

WHAT IS STILL MISSING, and it is now a short step rather than a vague one
-------------------------------------------------------------------------
The Planck SPECTRAL DISTRIBUTION is not assembled anywhere in the book.  The pieces
are all present: bosonic statistics for integer-spin excitations, by CAR-composition
order (\S sec:spin-statistics-structural); masslessness, hence the linear dispersion;
and the two polarisations of (2).  Occupation number 1/(e^{hbar w / k_B T} - 1) times
those gives Planck, and integrating it returns exactly the a of (5), which is checked
below.  What is absent is the derivation written down, not an ingredient.

SUPERSEDED IN PART, by comp 161 (the assembly), and CORRECTED HERE.
-------------------------------------------------------------------
This docstring originally ended by saying that "the photoelectric effect and the
quantisation of the field into countable quanta are absent outright".  Half of that
was wrong.  The quantisation into countable quanta is NOT absent: \S sec:projected-qft
constructs the Fock space by canonical quantisation of the projected Lorentzian
action, and the occupation number is the eigenvalue of N = a^dag a, an integer.  It is
INHERITED rather than derived, which is a different and much weaker complaint, and the
book says so in those words.  The claim was written here without checking that
chapter, and it propagated to the body, the Appendix B row and the website before
comp 161 caught it.

The PHOTOELECTRIC EFFECT is genuinely absent, and remains so after comp 161: Planck's
law is a statement about the equilibrium field, whereas the photoelectric effect is a
statement about matter absorbing one quantum and needs the interaction vertex.

A "pass": (2) 2 states for null and 3 for timelike, over 200 momenta each;
(4) T proportional to rho^(1/4) symbolically; (5) the g = 2 radiation constant within
1e-4 of the measured value; and the Planck integral reproducing it.
"""

from __future__ import annotations
import numpy as np
import sympy as sp
from scipy.integrate import quad

ETA = np.diag([1.0, -1.0, -1.0, -1.0])

HBAR = 1.054571817e-34
C_L = 2.99792458e8
K_B = 1.380649e-23
A_MEASURED = 7.565733e-16          # radiation constant, J m^-3 K^-4


def polarisation_count(rng, null: bool, trials: int = 200):
    """dim{eps : k.eps = 0} minus the gauge direction, when the gauge direction is
    itself transverse (which happens exactly when k is null)."""
    dims = set()
    for _ in range(trials):
        d = rng.normal(size=3)
        d /= np.linalg.norm(d)
        k = np.array([1.0, *d]) if null else np.array([1.0, *(d / 3.0)])
        transverse = 4 - np.linalg.matrix_rank((ETA @ k).reshape(1, 4))
        gauge = 1 if abs(float(k @ ETA @ k)) < 1e-12 else 0
        dims.add(transverse - gauge)
    return dims


def stefan_boltzmann_symbolic():
    """p = rho/3 plus the first law forces T ~ rho^(1/4)."""
    rho = sp.Symbol('rho', positive=True)
    V = sp.Symbol('V', positive=True)
    T = sp.Function('T', positive=True)(rho)
    # dS = (V/T) drho + (4 rho)/(3T) dV ; equality of mixed second partials
    resid = sp.simplify(sp.diff(V / T, V) - sp.diff(4 * rho / (3 * T), rho))
    sol = sp.dsolve(sp.numer(sp.together(resid)), T)
    return sol


def planck_integral(g=2):
    """Integrate the Planck distribution to get the radiation constant."""
    def f(x):
        if x > 700.0:                                # expm1 overflows; integrand is nil
            return 0.0
        return x ** 3 / np.expm1(x) if x > 0 else 0.0
    I, _ = quad(f, 0, np.inf, limit=200)             # = pi^4/15
    return g * K_B ** 4 * I / (2 * np.pi ** 2 * HBAR ** 3 * C_L ** 3)


def main():
    print("=" * 100)
    print("  Computation 160, the photon: what tau_null fixes, and what is still missing")
    print("=" * 100)
    print()
    rng = np.random.default_rng(20260927)

    print("  (2) EXACTLY TWO POLARISATIONS, and the discriminant is the null condition")
    dn = polarisation_count(rng, null=True)
    dt = polarisation_count(rng, null=False)
    print(f"      null momentum     (k.k = 0, the photon) : physical states = {sorted(dn)}")
    print(f"      timelike momentum (k.k > 0, massive)    : physical states = {sorted(dt)}")
    print(f"      the gauge direction is transverse to itself only when k.k = 0, which is")
    print(f"      the condition eq:metric puts at T = tau_null")
    check2 = dn == {2} and dt == {3}
    print(f"      -> two for the photon, three for a massive vector: {check2}")
    print()

    print("  (4) STEFAN-BOLTZMANN FROM w = 1/3 AND THE FIRST LAW, no quantum input")
    sol = stefan_boltzmann_symbolic()
    print(f"      integrability of dS with p = rho/3 gives  {sol}")
    exponent = sp.Rational(1, 4)
    # A substring test on a printer's output would match 11/4 or 1/40 and breaks on
    # any change of printing form.  Read the exponent of rho out of the solution
    # symbolically instead, and compare it with the value this line constructs.
    rho_sym = sp.Symbol('rho', positive=True)
    powers = [e.exp for e in sol.rhs.atoms(sp.Pow) if e.base == rho_sym]
    check4 = len(powers) == 1 and sp.simplify(powers[0] - exponent) == 0
    print(f"      so T ~ rho^{exponent}, hence rho ~ T^4: {check4}")
    print()

    print("  (5) AND THE TWO OF (2) FIXES THE CONSTANT")
    a2, a3 = planck_integral(2), planck_integral(3)
    print(f"      g = 2 polarisations : a = {a2:.6e} J m^-3 K^-4")
    print(f"      measured            : a = {A_MEASURED:.6e}")
    print(f"      relative difference : {abs(a2-A_MEASURED)/A_MEASURED:.2e}")
    print(f"      g = 3 would give    : a = {a3:.6e}, wrong by {abs(a3-A_MEASURED)/A_MEASURED:.0%}")
    check5 = abs(a2 - A_MEASURED) / A_MEASURED < 1e-4
    print(f"      -> the null condition fixes the coefficient, not only the T^4 scaling: {check5}")
    print()

    print("  SCOPE: what is missing is the derivation, not an ingredient")
    print("      The Planck spectral distribution is nowhere in the book.  Its pieces are:")
    print("      bosonic statistics for integer-spin excitations by CAR-composition order,")
    print("      masslessness hence linear dispersion, and the two polarisations above.")
    print("      comp 161 assembles it and returns exactly the constant of (5).")
    print("      CORRECTION: this computation first said the quantisation into countable")
    print("      quanta was absent outright.  It is not; \\S sec:projected-qft has the Fock")
    print("      space by canonical quantisation, inherited rather than derived.  The")
    print("      photoelectric effect IS absent, and stays absent after comp 161.")
    print()

    print("=" * 100)
    print("  ASSESSMENT")
    print("=" * 100)
    print(f"  (2) two polarisations for null, three for timelike            : {check2}")
    print(f"  (4) rho ~ T^4 from w = 1/3 and the first law alone            : {check4}")
    print(f"  (5) g = 2 reproduces the measured radiation constant          : {check5}")
    print()
    ok = check2 and check4 and check5
    print("  RESULT: " + ("ONE OBJECT, tau_null, fixes the light cone, the two polarisations, "
                          "the equation of state and the Stefan-Boltzmann law including its "
                          "coefficient; the Planck spectral distribution is assembled in comp 161"
                          if ok else "MISMATCH -- revise"))


if __name__ == "__main__":
    main()
