Black holes, rendered

"""A Schwarzschild black hole with a thin, opaque, glowing dust disc, seen from a camera near the disc's plane
so the far side of the disc is lensed over and under the shadow (the owner's ask, 2026-09-08).

    blackhole.py <outdir> [frames] [width] [height]

Geometric units, M = 1 (the horizon at r = 2, the photon sphere at 3, the innermost stable orbit at 6). The
camera sits at distance D = 40 at inclination 78 degrees from the disc's axis. For every pixel the null
geodesic is integrated once in its own orbital plane with the Binet equation u'' = -u + 3 u^2 (u = 1 / r),
fourth-order Runge-Kutta, until it falls through the horizon (black), escapes (a faint lensed star field
looked up by its exit direction) or crosses the disc's plane between r_in = 6 and r_out = 18 (a hit: the
disc's radius and azimuth there, and the redshift factor g = sqrt(1 - 3/r) / (1 + Omega b_z) for a Keplerian
emitter with Omega = r^-1.5 and b_z the photon's angular momentum about the disc's axis, which carries the
gravitational and Doppler shifts together). Only the first crossing counts: the disc is opaque.

Then every frame is a lookup: emission = g^4 * (1 - sqrt(6 / r)) r^-3 * dust(r, phi - omega(r) t), the dust a
multi-octave noise in (log r, phi) that winds up under differential rotation, and a warm colour ramp by the
observed temperature. Output: rgb24 raw frames for ffmpeg (<outdir>/frames.rgb), plus one PNG-ready preview.
The tracing takes a few minutes at 1280x720; the frames are fast.
"""
import math
import pathlib
import sys
import time

import numpy as np

out = pathlib.Path(sys.argv[1]); out.mkdir(parents=True, exist_ok=True)
NF = int(sys.argv[2]) if len(sys.argv) > 2 else 240
W = int(sys.argv[3]) if len(sys.argv) > 3 else 1280
H = int(sys.argv[4]) if len(sys.argv) > 4 else 720
FPS = 24.0

D = 40.0                      # camera distance
INC = math.radians(78.0)      # inclination from the disc's axis
FOV = math.radians(56.0)      # horizontal field of view
R_IN, R_OUT = 6.0, 18.0
DPHI = 0.004                  # integration step in the orbital angle
PHI_MAX = 3.5 * math.pi       # rays that loop more than this are lost (they are the photon ring's tail)

t0 = time.time()
# ---- camera rays -------------------------------------------------------------------------------------------
c_hat = np.array([math.sin(INC), 0.0, math.cos(INC)])            # the camera's position direction
fwd = -c_hat
right = np.cross(fwd, np.array([0.0, 0.0, 1.0])); right /= np.linalg.norm(right)
up = np.cross(right, fwd)
tanx = math.tan(FOV / 2); tany = tanx * H / W
px = (np.arange(W) + 0.5) / W * 2 - 1
py = 1 - (np.arange(H) + 0.5) / H * 2
PX, PY = np.meshgrid(px, py)
d = fwd[None, None, :] + tanx * PX[..., None] * right[None, None, :] + tany * PY[..., None] * up[None, None, :]
d = d.reshape(-1, 3); d /= np.linalg.norm(d, axis=1, keepdims=True)
N = d.shape[0]
e1 = np.repeat(c_hat[None, :], N, axis=0)                        # the orbital plane: e1 toward the camera
d_r = d @ c_hat
e2 = d - d_r[:, None] * e1[:, :]
d_t = np.linalg.norm(e2, axis=1)
e2 /= np.maximum(d_t, 1e-9)[:, None]
b_z = D * np.cross(c_hat, d)[:, 2]                               # the photon's L_z / E (sign chooses the disc's spin)

u = np.full(N, 1.0 / D, np.float64)
du = -u * d_r / np.maximum(d_t, 1e-9)                            # du/dphi at the camera (inward: positive)
phi = np.zeros(N)
z_prev = D * c_hat[2] * np.ones(N)
state = np.zeros(N, np.int8)                                     # 0 flying, 1 hit the disc, 2 fell in, 3 escaped
hit_r = np.zeros(N); hit_az = np.zeros(N); esc_dir = np.zeros((N, 3))
e1z, e2z = c_hat[2], e2[:, 2]

def f(u, du):
    return du, -u + 3.0 * u * u

active = np.arange(N)
step = 0
while active.size and step * DPHI < PHI_MAX:
    ua, dua = u[active], du[active]
    k1u, k1d = f(ua, dua)
    k2u, k2d = f(ua + 0.5 * DPHI * k1u, dua + 0.5 * DPHI * k1d)
    k3u, k3d = f(ua + 0.5 * DPHI * k2u, dua + 0.5 * DPHI * k2d)
    k4u, k4d = f(ua + DPHI * k3u, dua + DPHI * k3d)
    un = ua + DPHI / 6 * (k1u + 2 * k2u + 2 * k3u + k4u)
    dn = dua + DPHI / 6 * (k1d + 2 * k2d + 2 * k3d + k4d)
    phin = phi[active] + DPHI
    r = 1.0 / np.maximum(un, 1e-9)
    z = r * (math.cos(0) * 0 + np.cos(phin) * e1z + np.sin(phin) * e2z[active])
    fell = un >= 0.5
    escaped = (un < 1.0 / (2.5 * D)) & (dn < 0)
    crossed = (np.sign(z) != np.sign(z_prev[active])) & (z_prev[active] != 0) & ~fell
    # the crossing point, interpolated in phi
    frac = np.where(crossed, np.abs(z_prev[active]) / np.maximum(np.abs(z_prev[active]) + np.abs(z), 1e-9), 0.0)
    phic = phi[active] + frac * DPHI
    uc = u[active] + frac * (un - u[active])
    rc = 1.0 / np.maximum(uc, 1e-9)
    hit = crossed & (rc >= R_IN) & (rc <= R_OUT)
    # azimuth in the disc plane at the crossing
    pos = rc[:, None] * (np.cos(phic)[:, None] * e1[active] + np.sin(phic)[:, None] * e2[active])
    az = np.arctan2(pos[:, 1], pos[:, 0])
    idx = active
    state[idx[hit]] = 1; hit_r[idx[hit]] = rc[hit]; hit_az[idx[hit]] = az[hit]
    state[idx[fell & ~hit]] = 2
    esc = escaped & ~hit & ~fell
    state[idx[esc]] = 3
    esc_dir[idx[esc]] = (np.cos(phin)[:, None] * e1[active] + np.sin(phin)[:, None] * e2[active])[esc]
    u[idx] = un; du[idx] = dn; phi[idx] = phin; z_prev[idx] = z
    active = idx[~(hit | fell | esc)]
    step += 1
    if step % 400 == 0:
        print(f"  step {step}: {active.size} rays flying, {int((state == 1).sum())} hits, {int((state == 2).sum())} fell in, {int((state == 3).sum())} escaped ({time.time() - t0:.0f} s)")
state[active] = 2                                                # what is still circling is the photon ring's tail: black
print(f"traced {N} rays in {time.time() - t0:.0f} s: {int((state == 1).sum())} hit the disc, {int((state == 2).sum())} fell in, {int((state == 3).sum())} escaped")

# ---- the redshift, the geometry maps -----------------------------------------------------------------------
hitm = state == 1
omega = np.where(hitm, hit_r ** -1.5, 0.0)
g = np.where(hitm, np.sqrt(np.clip(1.0 - 3.0 / np.maximum(hit_r, 3.01), 0, 1)) / (1.0 + omega * b_z), 0.0)
emis = np.where(hitm, (1.0 - np.sqrt(R_IN / np.maximum(hit_r, R_IN))) * np.maximum(hit_r, R_IN) ** -3.0, 0.0)
emis /= emis[hitm].max() if hitm.any() else 1.0

# ---- the star field for escaped rays: sparse points by a hash of the exit direction -------------------------
def stars(dirs):
    v = np.floor((dirs + 1.0) * 400.0).astype(np.int64)
    h = (v[:, 0] * 73856093 ^ v[:, 1] * 19349663 ^ v[:, 2] * 83492791) & 0xFFFFFF
    s = (h < 0xFFFFFF * 0.004).astype(np.float64) * (0.3 + 0.7 * ((h >> 8) & 0xFF) / 255.0)
    return s
sky = np.zeros(N); escm = state == 3
sky[escm] = stars(esc_dir[escm]) * 0.9

# ---- the dust: multi-octave noise in (log r, phi), periodic in phi ------------------------------------------
rng = np.random.default_rng(3)
NP, NR = 512, 128
def band_noise(nr, nph, k_lo, k_hi):
    fr = np.fft.fftfreq(nr)[:, None]; fp = np.fft.fftfreq(nph)[None, :]
    k = np.hypot(fr * nr / 8.0, fp * nph / 40.0)
    band = (k >= k_lo) & (k <= k_hi)
    spec = np.where(band, np.exp(1j * rng.uniform(0, 2 * np.pi, (nr, nph))), 0.0)
    x = np.fft.ifft2(spec).real
    return (x - x.mean()) / (x.std() + 1e-9)
dust = 0.6 * band_noise(NR, NP, 0.5, 2.0) + 0.3 * band_noise(NR, NP, 2.0, 6.0) + 0.15 * band_noise(NR, NP, 6.0, 16.0)
dust = np.clip(0.5 + 0.45 * dust, 0.03, 1.4)

def ramp(x):
    """a warm ramp: black -> deep red -> orange -> yellow -> white, x in [0, ~1.5]"""
    x = np.clip(x, 0, 1.6)
    r = np.clip(x * 1.6, 0, 1)
    gch = np.clip(x * 1.15 - 0.25, 0, 1) ** 1.3
    b = np.clip(x * 0.9 - 0.55, 0, 1) ** 1.6
    return np.stack([r, gch, b], -1)

lr = np.where(hitm, (np.log(np.maximum(hit_r, R_IN)) - math.log(R_IN)) / (math.log(R_OUT) - math.log(R_IN)), 0.0)
ri = np.clip((lr * (NR - 1)).astype(int), 0, NR - 1)
omega_vis = np.where(hitm, 2 * np.pi / 4.0 * (R_IN / np.maximum(hit_r, R_IN)) ** 1.5, 0.0)   # one turn in 4 s at r_in
sky_rgb = np.stack([sky, sky, sky * 1.1], -1)

frames = np.zeros((NF, H * W, 3), np.uint8)
for n in range(NF):
    t = n / FPS
    ph = np.mod(hit_az - omega_vis * t, 2 * np.pi)
    pi_ = np.clip((ph / (2 * np.pi) * NP).astype(int), 0, NP - 1)
    tex = dust[ri, pi_]
    I = emis * g ** 4 * tex * 6.0
    temp = np.where(hitm, g * (R_IN / np.maximum(hit_r, R_IN)) ** 0.75, 0.0)
    col = ramp(temp * 1.15) * (I / (1.0 + I))[:, None] * 1.5 + sky_rgb
    frames[n] = (np.clip(col, 0, 1) ** (1 / 2.2) * 255).astype(np.uint8)
    if n % 60 == 0:
        print(f"  frame {n} ({time.time() - t0:.0f} s)")
frames.tofile(out / "frames.rgb")
print(f"{NF} frames of {W}x{H} -> {out / 'frames.rgb'} ({time.time() - t0:.0f} s)")

"""A black hole that is not Schwarzschild: a thin glowing dust disc around a Kerr black hole, or around a
Johannsen-Psaltis deformation of one that solves no vacuum field equation, rendered by one metric-agnostic
null-geodesic tracer (2026-09-08).

    geodesic_disc.py <metric> <outdir> [frames] [width] [height] [bands] [spin] [eps3]
        metric  schwarzschild | kerr | jp
        frames  default 120 at 24 fps
        width, height  default 960 x 540 (3840 x 2160 takes about two hours of CPU with 4 bands)
        bands   horizontal bands the frame is traced in, default 1 (use 4 for 4K: memory); each band is
                checkpointed in <outdir>, so an interrupted run resumes where it stopped
        spin    a / M for kerr and jp, default 0.9
        eps3    the deformation for jp, default 3.0 (0 is Kerr again)

THE TRACER. Every pixel's photon is integrated backwards from the camera with Hamilton's equations for a null
geodesic, H = 1/2 g^{mu nu} p_mu p_nu = 0, in Boyer-Lindquist-like coordinates (t, r, theta, phi):
    dx^mu / dlambda = g^{mu nu} p_nu,      dp_r / dlambda = -dH/dr,      dp_theta / dlambda = -dH/dtheta,
with p_t = -1 (energy 1) and p_phi conserved. The metric enters ONLY as its covariant components g_tt, g_tphi,
g_rr, g_thetatheta, g_phiphi as functions of (r, theta): the inverse is taken numerically (a 2x2 block and two
reciprocals) and dH/dr, dH/dtheta by central differences, so any stationary axisymmetric metric can be dropped
into metric_cov() without deriving anything. Fourth-order Runge-Kutta, the step shrinking toward the horizon
and toward the coordinate poles (without the latter a seam appears at the top of the shadow, where rays pass
over the pole). A ray ends when it falls through the horizon (black), escapes beyond 55 M (a star field by its
exit direction), or crosses the disc's plane theta = pi/2 between the innermost stable circular orbit and
18 M (a hit). The camera is a zero-angular-momentum observer at 40 M, 78 degrees off the disc's axis.

THE METRICS (geometric units, M = 1).
    schwarzschild   a = 0: the control (horizon 2, innermost stable orbit 6).
    kerr            frame dragging, the D-shaped shadow flattened on the prograde side, the disc reaching in
                    (a = 0.9: horizon 1.436, innermost stable orbit 2.321).
    jp              Johannsen & Psaltis (2011): Kerr with the deformation h = eps3 M^3 r / Sigma^2 in
                    g_tt, g_tphi, g_rr and g_phiphi. NOT a solution of the vacuum Einstein equations: a
                    parametrised deviation of the kind used to test whether real black holes are Kerr. The
                    image is what integrating this metric's geodesics gives, nothing else (a = 0.9, eps3 = 3:
                    innermost stable orbit 1.465, just outside the horizon).

THE DISC. Thin, opaque, on prograde circular orbits. The orbital frequency Omega(r), u^t and the innermost
stable orbit are found numerically from the metric in the equatorial plane (Omega from the radial derivatives
of g_tt, g_tphi, g_phiphi; the innermost stable orbit at the minimum of the orbit's energy), so a new metric
brings its own disc. The redshift factor of a hit is Cunningham's g = 1 / (u^t (1 - Omega p_phi)) with p_phi
the photon's conserved angular momentum: gravitational and Doppler shifts together. The emission is
g^4 (1 - sqrt(r_in / r)) r^-3 times a multi-octave dust pattern in (log r, phi) that winds up under the
metric's own Omega(r) (one turn in 4 s at the inner edge), coloured by the observed temperature on a warm
ramp. Output: <outdir>/frames.rgb, raw rgb24 frames for ffmpeg, e.g.
    ffmpeg -f rawvideo -pix_fmt rgb24 -s 960x540 -r 24 -i frames.rgb -vf format=yuv420p -c:v libx264 out.mp4
and <outdir>/summary.txt with the horizon, the innermost stable orbit and the ray counts.

This is not part of the interpolation tests. It was written for the owner's curiosity and stays because the
tracer is general and small.
"""
import math
import pathlib
import sys
import time

import numpy as np

metric_name = sys.argv[1]
out = pathlib.Path(sys.argv[2]); out.mkdir(parents=True, exist_ok=True)
NF = int(sys.argv[3]) if len(sys.argv) > 3 else 120
W = int(sys.argv[4]) if len(sys.argv) > 4 else 960
H = int(sys.argv[5]) if len(sys.argv) > 5 else 540
BANDS = int(sys.argv[6]) if len(sys.argv) > 6 else 1
SPIN = float(sys.argv[7]) if len(sys.argv) > 7 else 0.9
EPS3_ARG = float(sys.argv[8]) if len(sys.argv) > 8 else 3.0
FPS = 24.0
if metric_name not in ("schwarzschild", "kerr", "jp"):
    sys.exit("metric: schwarzschild | kerr | jp")
A = 0.0 if metric_name == "schwarzschild" else SPIN
EPS3 = EPS3_ARG if metric_name == "jp" else 0.0
R_OBS, TH_OBS = 40.0, math.radians(78.0)
FOV = math.radians(56.0)
R_OUT = 18.0
R_ESC = 55.0
t0 = time.time()


def metric_cov(r, th):
    """Covariant components g_tt, g_tp, g_rr, g_hh, g_pp of the Kerr metric with the Johannsen-Psaltis
    deformation h = EPS3 r / Sigma^2 (EPS3 = 0 is Kerr; A = 0 too is Schwarzschild). M = 1."""
    s2 = np.maximum(np.sin(th) ** 2, 1e-12)
    c2 = 1.0 - s2
    Sig = r * r + A * A * c2
    Del = r * r - 2.0 * r + A * A
    h = EPS3 * r / (Sig * Sig)
    g_tt = -(1.0 + h) * (1.0 - 2.0 * r / Sig)
    g_tp = -2.0 * A * r * s2 * (1.0 + h) / Sig
    g_pp = s2 * (r * r + A * A + 2.0 * A * A * r * s2 / Sig) + h * A * A * s2 * s2 * (Sig + 2.0 * r) / Sig
    g_rr = Sig * (1.0 + h) / np.maximum(Del + A * A * h * s2, 1e-9)
    g_hh = Sig
    return g_tt, g_tp, g_rr, g_hh, g_pp


def metric_inv(r, th):
    g_tt, g_tp, g_rr, g_hh, g_pp = metric_cov(r, th)
    det = g_tt * g_pp - g_tp * g_tp
    det = np.where(np.abs(det) < 1e-12, -1e-12, det)
    return g_pp / det, -g_tp / det, 1.0 / g_rr, 1.0 / g_hh, g_tt / det      # g^tt, g^tp, g^rr, g^hh, g^pp


def hamiltonian(r, th, p_t, p_r, p_h, p_p):
    gtt, gtp, grr, ghh, gpp = metric_inv(r, th)
    return 0.5 * (gtt * p_t * p_t + 2.0 * gtp * p_t * p_p + gpp * p_p * p_p + grr * p_r * p_r + ghh * p_h * p_h)


def rhs(r, th, p_r, p_h, p_t, p_p):
    gtt, gtp, grr, ghh, gpp = metric_inv(r, th)
    dr = grr * p_r
    dth = ghh * p_h
    dph = gtp * p_t + gpp * p_p
    hr = 1e-4 * np.maximum(r, 1.0)
    dHdr = (hamiltonian(r + hr, th, p_t, p_r, p_h, p_p) - hamiltonian(r - hr, th, p_t, p_r, p_h, p_p)) / (2.0 * hr)
    hh = 1e-4
    dHdh = (hamiltonian(r, th + hh, p_t, p_r, p_h, p_p) - hamiltonian(r, th - hh, p_t, p_r, p_h, p_p)) / (2.0 * hh)
    return dr, dth, dph, -dHdr, -dHdh


# ---- the horizon and the disc's orbits from the metric ------------------------------------------------------
r_h = 1.0 + math.sqrt(max(1.0 - A * A, 0.0))
rr = np.linspace(r_h * 1.02, 30.0, 6000)
th_eq = np.full_like(rr, math.pi / 2)
dr_ = 1e-4
def eq(r):
    return metric_cov(r, np.full_like(r, math.pi / 2))
gtt0, gtp0, _, _, gpp0 = eq(rr)
gtt1, gtp1, _, _, gpp1 = eq(rr + dr_); gtt2, gtp2, _, _, gpp2 = eq(rr - dr_)
dgtt = (gtt1 - gtt2) / (2 * dr_); dgtp = (gtp1 - gtp2) / (2 * dr_); dgpp = (gpp1 - gpp2) / (2 * dr_)
disc_ = dgtp * dgtp - dgtt * dgpp
Om = (-dgtp + np.sqrt(np.maximum(disc_, 0.0))) / dgpp                       # prograde
norm = -(gtt0 + 2.0 * gtp0 * Om + gpp0 * Om * Om)
ok = (norm > 0) & (disc_ > 0)
E_orb = np.where(ok, -(gtt0 + gtp0 * Om) / np.sqrt(np.maximum(norm, 1e-12)), np.inf)
i_isco = int(np.argmin(E_orb))
R_IN = float(rr[i_isco])
def orbit(r):
    """Omega and u^t of the prograde circular orbit at equatorial r, by interpolation on the grid."""
    Omg = np.interp(r, rr, Om)
    ut = np.interp(r, rr, 1.0 / np.sqrt(np.maximum(norm, 1e-12)))
    return Omg, ut
print(f"{metric_name}: a = {A}, eps3 = {EPS3}, horizon r = {r_h:.3f}, ISCO r = {R_IN:.3f} (Omega there {float(Om[i_isco]):.4f})")

# ---- the camera: a zero-angular-momentum observer's tetrad at (R_OBS, TH_OBS, phi = 0) ---------------------
g_tt, g_tp, g_rr, g_hh, g_pp = [float(x[0]) for x in metric_cov(np.array([R_OBS]), np.array([TH_OBS]))]
omega_z = -g_tp / g_pp
alpha = math.sqrt(-(g_tt - g_tp * g_tp / g_pp))
tanx = math.tan(FOV / 2); tany = tanx * H / W
px = (np.arange(W) + 0.5) / W * 2 - 1
py = 1 - (np.arange(H) + 0.5) / H * 2


def trace_band(y0, y1):
    """Trace the rays of rows y0..y1-1; return state, hit_r, hit_az, hit_pp, esc_dir for them."""
    PX, PY = np.meshgrid(px, py[y0:y1])
    PX = PX.ravel(); PY = PY.ravel()
    N = PX.size
    n_r = -np.ones(N); n_h = -tany * PY; n_p = -tanx * PX
    k = 1.0 / np.sqrt(n_r * n_r + n_h * n_h + n_p * n_p)
    n_r *= k; n_h *= k; n_p *= k
    pu_t = 1.0 / alpha * np.ones(N)
    pu_r = n_r / math.sqrt(g_rr)
    pu_h = n_h / math.sqrt(g_hh)
    pu_p = omega_z / alpha + n_p / math.sqrt(g_pp)
    p_t = g_tt * pu_t + g_tp * pu_p
    p_p = g_tp * pu_t + g_pp * pu_p
    p_r = g_rr * pu_r
    p_h = g_hh * pu_h
    E = -p_t
    p_t = p_t / E; p_p = p_p / E; p_r = p_r / E; p_h = p_h / E
    r = np.full(N, R_OBS); th = np.full(N, TH_OBS); ph = np.zeros(N)
    state = np.zeros(N, np.int8)
    hit_r = np.zeros(N); hit_az = np.zeros(N); hit_pp = np.zeros(N); esc_dir = np.zeros((N, 3))
    active = np.arange(N)
    step = 0
    while active.size and step < MAX_STEPS:
        ra, ta, pa, pra, pha = r[active], th[active], ph[active], p_r[active], p_h[active]
        pt, pp = p_t[active], p_p[active]
        dl = np.clip(0.04 * (ra - r_h), 0.004, 0.5) * np.clip(np.sin(ta) / 0.25, 0.04, 1.0)
        k1 = rhs(ra, ta, pra, pha, pt, pp)
        k2 = rhs(ra + 0.5 * dl * k1[0], ta + 0.5 * dl * k1[1], pra + 0.5 * dl * k1[3], pha + 0.5 * dl * k1[4], pt, pp)
        k3 = rhs(ra + 0.5 * dl * k2[0], ta + 0.5 * dl * k2[1], pra + 0.5 * dl * k2[3], pha + 0.5 * dl * k2[4], pt, pp)
        k4 = rhs(ra + dl * k3[0], ta + dl * k3[1], pra + dl * k3[3], pha + dl * k3[4], pt, pp)
        rn = ra + dl / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0])
        tn = ta + dl / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1])
        pn = pa + dl / 6 * (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2])
        prn = pra + dl / 6 * (k1[3] + 2 * k2[3] + 2 * k3[3] + k4[3])
        phn = pha + dl / 6 * (k1[4] + 2 * k2[4] + 2 * k3[4] + k4[4])
        fell = rn <= r_h * 1.02
        escaped = (rn > R_ESC) & (rn > ra)
        zc = (np.cos(ta) * np.cos(tn) < 0)
        frac = np.where(zc, np.abs(np.cos(ta)) / np.maximum(np.abs(np.cos(ta)) + np.abs(np.cos(tn)), 1e-12), 0.0)
        rc = ra + frac * (rn - ra); pc = pa + frac * (pn - pa)
        hit = zc & (rc >= R_IN) & (rc <= R_OUT) & ~fell
        idx = active
        state[idx[hit]] = 1; hit_r[idx[hit]] = rc[hit]; hit_az[idx[hit]] = pc[hit]; hit_pp[idx[hit]] = pp[hit]
        state[idx[fell & ~hit]] = 2
        esc = escaped & ~hit & ~fell
        state[idx[esc]] = 3
        esc_dir[idx[esc]] = np.stack([np.sin(tn) * np.cos(pn), np.sin(tn) * np.sin(pn), np.cos(tn)], -1)[esc]
        r[idx] = rn; th[idx] = np.clip(tn, 1e-4, math.pi - 1e-4); ph[idx] = pn; p_r[idx] = prn; p_h[idx] = phn
        active = idx[~(hit | fell | esc)]
        step += 1
    state[active] = 2
    return state, hit_r, hit_az, hit_pp, esc_dir


MAX_STEPS = 40000
parts = []
for b in range(BANDS):
    y0, y1 = b * H // BANDS, (b + 1) * H // BANDS
    ck = out / f"band_{b}_of_{BANDS}_{W}x{H}.npz"          # each band is checkpointed: a restart skips the bands done
    if ck.exists():
        z = np.load(ck)
        parts.append((z["state"], z["hit_r"], z["hit_az"], z["hit_pp"], z["esc_dir"]))
        print(f"  band {b + 1}/{BANDS} (rows {y0}-{y1 - 1}): loaded from {ck.name}", flush=True)
        continue
    parts.append(trace_band(y0, y1))
    st, hr_, ha_, hp_, ed_ = parts[-1]
    np.savez(ck, state=st, hit_r=hr_, hit_az=ha_, hit_pp=hp_, esc_dir=ed_)
    print(f"  band {b + 1}/{BANDS} (rows {y0}-{y1 - 1}): {int((st == 1).sum())} hit, {int((st == 2).sum())} fell, {int((st == 3).sum())} escaped ({time.time() - t0:.0f} s)", flush=True)
state = np.concatenate([p[0] for p in parts]); hit_r = np.concatenate([p[1] for p in parts]); hit_az = np.concatenate([p[2] for p in parts])
hit_pp = np.concatenate([p[3] for p in parts]); esc_dir = np.concatenate([p[4] for p in parts]); del parts
N = state.size
print(f"traced {N} rays in {time.time() - t0:.0f} s: {int((state == 1).sum())} hit, {int((state == 2).sum())} fell, {int((state == 3).sum())} escaped", flush=True)

# ---- the redshift and the frames ---------------------------------------------------------------------------
hitm = state == 1
Omg, ut = orbit(np.where(hitm, hit_r, R_IN))
g = np.where(hitm, 1.0 / (ut * (1.0 - Omg * hit_pp)), 0.0)
g = np.clip(g, 0.0, 3.0)
emis = np.where(hitm, (1.0 - np.sqrt(R_IN / np.maximum(hit_r, R_IN))) * np.maximum(hit_r, R_IN) ** -3.0, 0.0)
emis /= emis[hitm].max() if hitm.any() else 1.0

def stars(dirs):
    v = np.floor((dirs + 1.0) * (400.0 * W / 960.0)).astype(np.int64)
    hsh = (v[:, 0] * 73856093 ^ v[:, 1] * 19349663 ^ v[:, 2] * 83492791) & 0xFFFFFF
    return (hsh < 0xFFFFFF * 0.004 * (960.0 / W) ** 0.5).astype(np.float64) * (0.3 + 0.7 * ((hsh >> 8) & 0xFF) / 255.0)
sky = np.zeros(N); escm = state == 3
sky[escm] = stars(esc_dir[escm]) * 0.9

rng = np.random.default_rng(3)
NP = int(512 * W / 960); NR = int(128 * W / 960)
def band_noise(nr, nph, k_lo, k_hi):
    fr = np.fft.fftfreq(nr)[:, None]; fp = np.fft.fftfreq(nph)[None, :]
    kk = np.hypot(fr * nr / (8.0 * W / 960.0), fp * nph / (40.0 * W / 960.0))
    band = (kk >= k_lo) & (kk <= k_hi)
    spec = np.where(band, np.exp(1j * rng.uniform(0, 2 * np.pi, (nr, nph))), 0.0)
    x = np.fft.ifft2(spec).real
    return (x - x.mean()) / (x.std() + 1e-9)
dust = 0.6 * band_noise(NR, NP, 0.5, 2.0) + 0.3 * band_noise(NR, NP, 2.0, 6.0) + 0.15 * band_noise(NR, NP, 6.0, 16.0)
dust = np.clip(0.5 + 0.45 * dust, 0.03, 1.4)

def ramp(x):
    x = np.clip(x, 0, 1.6)
    return np.stack([np.clip(x * 1.6, 0, 1), np.clip(x * 1.15 - 0.25, 0, 1) ** 1.3, np.clip(x * 0.9 - 0.55, 0, 1) ** 1.6], -1)

lr = np.where(hitm, (np.log(np.maximum(hit_r, R_IN)) - math.log(R_IN)) / (math.log(R_OUT) - math.log(R_IN)), 0.0)
ri = np.clip((lr * (NR - 1)).astype(int), 0, NR - 1)
Om_in = float(np.interp(R_IN, rr, Om))
omega_vis = np.where(hitm, 2 * np.pi / 4.0 * Omg / Om_in, 0.0)        # the metric's own Omega(r), one turn in 4 s at the ISCO
temp_scale = np.where(hitm, g * (R_IN / np.maximum(hit_r, R_IN)) ** 0.75, 0.0)
sky_rgb = np.stack([sky, sky, sky * 1.1], -1)
base_col = ramp(temp_scale * 1.15) * 1.5
with open(out / "frames.rgb", "wb") as fh:
    for n in range(NF):
        t = n / FPS
        phn_ = np.mod(hit_az - omega_vis * t, 2 * np.pi)
        pi_ = np.clip((phn_ / (2 * np.pi) * NP).astype(int), 0, NP - 1)
        I = emis * g ** 4 * dust[ri, pi_] * 6.0
        col = base_col * (I / (1.0 + I))[:, None] + sky_rgb
        fh.write((np.clip(col, 0, 1) ** (1 / 2.2) * 255).astype(np.uint8).tobytes())
        if n % 24 == 0:
            print(f"  frame {n} ({time.time() - t0:.0f} s)", flush=True)
(out / "summary.txt").write_text(f"{metric_name} a={A} eps3={EPS3} horizon={r_h:.4f} isco={R_IN:.4f} hits={int(hitm.sum())} fell={int((state == 2).sum())} escaped={int(escm.sum())} g_max={float(g[hitm].max()) if hitm.any() else 0:.3f}\n")
print(f"{NF} frames of {W}x{H} -> {out / 'frames.rgb'} ({time.time() - t0:.0f} s)", flush=True)

# A black hole that is not Schwarzschild

One script, `geodesic_disc.py`, numpy only, with ffmpeg for the encode. It renders a
thin glowing dust disc around a black hole by tracing every pixel's light ray
backwards through the spacetime, and it does not care which spacetime: the metric
enters as five functions of radius and angle, and the tracer inverts and
differentiates them numerically. Three are built in.

- `schwarzschild`: the control. Horizon at 2 M, innermost stable orbit at 6 M.
- `kerr`: the spinning hole (default spin 0.9). Frame dragging, the shadow
  flattened into a D on its prograde side, the disc reaching in to 2.32 M.
- `jp`: the Johannsen-Psaltis deformation of that Kerr metric (default deviation
  3). It is not a solution of the vacuum Einstein equations at all. It is the
  kind of parametrised non-Kerr black hole astronomers test the no-hair theorem
  against, and its image is whatever integrating its light paths gives.

    ./geodesic_disc.py jp out 120 960 540
    ffmpeg -f rawvideo -pix_fmt rgb24 -s 960x540 -r 24 -i out/frames.rgb -vf format=yuv420p -c:v libx264 jp.mp4

A 960x540 frame traces in about six minutes on one CPU; 3840x2160 takes about two
hours and should be run with four bands (the sixth argument) so it fits in memory.
The disc's inner edge, its rotation and its redshift are all found from the metric,
so dropping a new metric into `metric_cov()` brings its own disc with it.

This lives beside the interpolation project rather than inside it. It was written in
an evening for the owner's curiosity, after the question "can you render a black
hole that is not the Schwarzschild metric", and it stays because the tracer is
general and small. There is no general solution to render: that is numerical
relativity on a supercomputer. This is two steps beyond Schwarzschild, the second of
them beyond any exact solution.

"""Hawking radiation of a Schwarzschild black hole in photons: the greybody factors, the spectrum, the power, and
the mode functions at the spectrum's peak (2026-09-08, the owner's question "can Hawking radiation be rendered").

Geometric units G = c = hbar = k_B = 1, M = 1: T_H = 1 / (8 pi) = 0.0398. For each multipole l the electromagnetic
(s = 1) Regge-Wheeler equation  d2psi/dr*2 + (w^2 - V) psi = 0,  V = (1 - 2/r) l (l + 1) / r^2,  r* = r + 2 ln(r/2 - 1)
is integrated inward from r = R_far with the pure outgoing wave psi = exp(i w r*) (this is the "up" mode: born at the
horizon, partly transmitted to infinity, partly reflected back in) down to r - 2 = 1e-7, where it is decomposed into
A exp(i w r*) + B exp(-i w r*). The transmission probability Gamma_l(w) = 1 / |A|^2 is the greybody factor (by
reciprocity the same as for a wave sent in from infinity), and psi / A is the mode itself, unit outgoing amplitude at the
horizon. The photon emission is then

    dN/dt dw = 2 * sum_l (2l + 1) Gamma_l(w) / (2 pi (exp(8 pi w) - 1)),   dE/dt dw = w dN/dt dw,

the 2 for the polarisations. Checks printed: Gamma -> 1 at high w for low l, the w^4 law at low w, and the total photon
power against Page 1976 (3.36e-5 hbar c^6 / G^2 M^2, 16.7 % of his 2.011e-4 total). Writes greybody.npz: the w grid,
Gamma[l, w], the spectra, and the l = 1, 2 up-modes on an r* grid at the spectrum's peak frequency.
"""
import math
import pathlib
import sys
import time

import numpy as np

here = pathlib.Path(__file__).resolve().parent
OUT = here / "greybody.npz"
LMAX = 8
R_FAR = 3000.0
X_HOR = 1e-7                  # r - 2 where the horizon decomposition is done
H = 0.05                      # step in r*
w = np.exp(np.linspace(math.log(0.01), math.log(1.5), 260))   # the frequency grid, w M
NW = w.size
T_H = 1.0 / (8 * math.pi)


def rstar_of_x(x):
    return 2.0 + x + 2.0 * np.log(x / 2.0)


def x_of_rstar(rs):
    """invert r* = 2 + x + 2 ln(x/2) by Newton (rs may be an array)"""
    rs = np.asarray(rs, dtype=np.float64)
    x = np.where(rs > 4, rs - 2.0, 2.0 * np.exp((rs - 2.0) / 2.0))
    for _ in range(60):
        f = 2.0 + x + 2.0 * np.log(x / 2.0) - rs
        x = np.maximum(x - f / (1.0 + 2.0 / x), 1e-300)
    return x


def integrate_up(l, ws, keep_below=None):
    """the up-mode for multipole l at frequencies ws (vector): inward RK4 from R_FAR. Returns Gamma (per w), and if
    keep_below is given (an r* value) the arrays (rstar, x, psi[w, i]) for r* < keep_below."""
    ll = l * (l + 1)
    x0 = R_FAR - 2.0
    rs0 = rstar_of_x(x0)
    rs_end = rstar_of_x(X_HOR)
    n = int(math.ceil((rs0 - rs_end) / H))
    h = -(rs0 - rs_end) / n
    rs = rs0
    x = np.float64(x0)
    psi = np.exp(1j * ws * rs0)
    dpsi = 1j * ws * psi
    keep_rs, keep_x, keep_psi = [], [], []

    def V_of_x(xx):
        return (xx / (xx + 2.0)) * ll / (xx + 2.0) ** 2

    def dx(xx):
        return xx / (xx + 2.0)

    for i in range(n):
        # x at the stage points (scalar, shared by all w)
        k1x = dx(x)
        x2 = x + 0.5 * h * k1x
        k2x = dx(x2)
        x3 = x + 0.5 * h * k2x
        k3x = dx(x3)
        x4 = x + h * k3x
        k4x = dx(x4)
        V1, V2, V4 = V_of_x(x), V_of_x(x2), V_of_x(x4)       # V at stage 2 and 3 are the same point in r*
        w2 = ws * ws
        k1p = dpsi;                     k1d = (V1 - w2) * psi
        p2 = psi + 0.5 * h * k1p;       d2 = dpsi + 0.5 * h * k1d
        k2p = d2;                       k2d = (V2 - w2) * p2
        p3 = psi + 0.5 * h * k2p;       d3 = dpsi + 0.5 * h * k2d
        k3p = d3;                       k3d = (V2 - w2) * p3
        p4 = psi + h * k3p;             d4 = dpsi + h * k3d
        k4p = d4;                       k4d = (V4 - w2) * p4
        psi = psi + h / 6 * (k1p + 2 * k2p + 2 * k3p + k4p)
        dpsi = dpsi + h / 6 * (k1d + 2 * k2d + 2 * k3d + k4d)
        x = x + h / 6 * (k1x + 2 * k2x + 2 * k3x + k4x)
        rs = rs + h
        if keep_below is not None and rs < keep_below:
            keep_rs.append(rs); keep_x.append(float(x)); keep_psi.append(psi.copy())
    # decompose at the horizon end: psi = A e^{i w r*} + B e^{-i w r*}
    A = 0.5 * (psi + dpsi / (1j * ws)) * np.exp(-1j * ws * rs)
    B = 0.5 * (psi - dpsi / (1j * ws)) * np.exp(1j * ws * rs)
    gamma = 1.0 / np.abs(A) ** 2
    flux_check = np.abs(A) ** 2 - np.abs(B) ** 2      # must be 1 (unit transmitted flux at infinity)
    if keep_below is None:
        return gamma, flux_check, A, B
    keep = (np.array(keep_rs[::-1]), np.array(keep_x[::-1]), np.array(keep_psi[::-1]) / A[None, :])
    return gamma, flux_check, A, B, keep


t0 = time.time()
Gamma = np.zeros((LMAX + 1, NW))
for l in range(1, LMAX + 1):
    g, fc, A, B = integrate_up(l, w)
    Gamma[l] = g
    print(f"l = {l}: Gamma at w = {w[0]:.3f}: {g[0]:.3e}, at w = {w[NW // 2]:.3f}: {g[NW // 2]:.4f}, at w = {w[-1]:.2f}: {g[-1]:.6f}; "
          f"flux conservation max |1 - (|A|^2 - |B|^2)| = {np.abs(1 - fc).max():.1e} ({time.time() - t0:.0f} s)", flush=True)

# checks: the w^4 law at low w for l = 1 (Gamma ~ c w^4), and the geometric-optics limit sum (2l+1) Gamma -> 27 w^2
lo = (w < 0.03)
slope = np.polyfit(np.log(w[lo]), np.log(Gamma[1][lo]), 1)[0]
S = np.sum((2 * np.arange(LMAX + 1) + 1)[:, None] * Gamma, axis=0)
print(f"low-w slope of Gamma_1: {slope:.3f} (theory 4)")
for wi in (0.3, 0.6, 1.0, 1.5):
    i = np.argmin(np.abs(w - wi))
    print(f"  sum (2l+1) Gamma at w = {w[i]:.2f}: {S[i]:.3f}  vs geometric 27 w^2 = {27 * w[i] ** 2:.3f}  (ratio {S[i] / (27 * w[i] ** 2):.3f})")

# the spectrum and the power
planck = 1.0 / np.expm1(8 * math.pi * w)
dNdw = 2.0 * S * planck / (2 * math.pi)
dEdw = w * dNdw
dEdw_geo = 2.0 * 27 * w ** 2 * planck / (2 * math.pi) * w
P = np.trapezoid(dEdw, w); P_geo = np.trapezoid(dEdw_geo, w); Ndot = np.trapezoid(dNdw, w)
P_geo_exact = 27 * math.pi ** 3 * T_H ** 4 / 15
ipk = int(np.argmax(dEdw)); ipk_n = int(np.argmax(dNdw))
print(f"photon power P M^2 = {P:.3e}  (Page 1976: 3.36e-5; geometric-optics blackbody {P_geo:.3e}, exact {P_geo_exact:.3e}; ratio P/P_geo = {P / P_geo_exact:.3f})")
print(f"photon rate  N M   = {Ndot:.3e} per unit time M;  mean photon energy {P / Ndot:.4f} = {P / Ndot / T_H:.2f} T_H")
print(f"energy spectrum peak at w M = {w[ipk]:.3f} = {w[ipk] / T_H:.2f} T_H  (blackbody peak 2.82 T_H = {2.82 * T_H:.3f});  wavelength 2 pi / w = {2 * math.pi / w[ipk]:.1f} M = {math.pi / w[ipk]:.1f} r_s")
print(f"number spectrum peak at w M = {w[ipk_n]:.3f}")
frac = {}
for l in range(1, 5):
    frac[l] = np.trapezoid(w * 2 * (2 * l + 1) * Gamma[l] * planck / (2 * math.pi), w) / P
print("share of the power by multipole: " + ", ".join(f"l={l}: {frac[l] * 100:.1f}%" for l in frac))

# the modes at the peak, l = 1 and 2, kept for r* < 90 (r < ~ 80)
w_pk = float(w[ipk])
modes = {}
for l in (1, 2):
    g, fc, A, B, (krs, kx, kpsi) = integrate_up(l, np.array([w_pk]), keep_below=90.0)
    modes[l] = (krs, kx, kpsi[:, 0], float(g[0]), complex(B[0] / A[0]))
    print(f"mode l = {l} at w = {w_pk:.3f}: Gamma = {g[0]:.4f}, |R|^2 = {abs(B[0] / A[0]) ** 2:.4f}, {krs.size} points kept, r* in [{krs[0]:.1f}, {krs[-1]:.1f}] ({time.time() - t0:.0f} s)")

np.savez(OUT, w=w, Gamma=Gamma, S=S, dNdw=dNdw, dEdw=dEdw, dEdw_geo=dEdw_geo, P=P, Ndot=Ndot, w_peak=w_pk,
         rs1=modes[1][0], x1=modes[1][1], psi1=modes[1][2], gamma1=modes[1][3],
         rs2=modes[2][0], x2=modes[2][1], psi2=modes[2][2], gamma2=modes[2][3])
print(f"wrote {OUT} ({time.time() - t0:.0f} s)")

"""Shared pieces for the Hawking renders (2026-09-08): colour from a spectrum, the Hawking spectrum from greybody.npz,
a 5x7 bitmap font for captions, and a few constants. numpy only.

Radiance bookkeeping (geometric units, M = 1, per unit frequency w per steradian, both polarisations):
    blackbody          B_w(T) = w^3 / (4 pi^3 (exp(w / T) - 1))
    Hawking, at infinity, along a ray from the horizon (the shadow's uniform radiance in ray optics):
                       I_w = [S(w) / (27 w^2)] * B_w(T_H),   S = sum_l (2l + 1) Gamma_l(w),   T_H = 1 / (8 pi)
    seen by a static observer at r (blueshift beta = 1 / sqrt(1 - 2/r), the same in every direction):
                       I_loc(w) = beta^3 I_inf(w / beta)      (I_w / w^3 is invariant along a ray)
A temperature in kelvin converts to geometric units through the hole's mass: T_geo = (T_K / T_H,K) / (8 pi), and a
frequency w (geometric) is a wavelength lambda = 2 pi c / (w c^3 / (G M)). The CIE 1931 colour-matching functions are
the multi-lobe Gaussian fits of Wyman, Sloan & Shirley (2013), then XYZ -> linear sRGB (D65).
"""
import math
import pathlib

import numpy as np

HBAR, C, G, KB = 1.0546e-34, 2.99792458e8, 6.674e-11, 1.3807e-23
T_H_GEO = 1.0 / (8 * math.pi)
MK_PER_TH = HBAR * C ** 3 / (8 * math.pi * G * KB)        # M * T_H in kg K: 1.227e23


def mass_for_TH(T_K):
    return MK_PER_TH / T_K


def rs_metres(M_kg):
    return 2 * G * M_kg / C ** 2


LAM = np.linspace(380e-9, 780e-9, 401)                     # the colour integral's wavelength grid


def _lobe(lam_nm, mu, s1, s2):
    s = np.where(lam_nm < mu, s1, s2)
    return np.exp(-0.5 * ((lam_nm - mu) / s) ** 2)


_ln = LAM * 1e9
XBAR = 1.056 * _lobe(_ln, 599.8, 37.9, 31.0) + 0.362 * _lobe(_ln, 442.0, 16.0, 26.7) - 0.065 * _lobe(_ln, 501.1, 20.4, 26.2)
YBAR = 0.821 * _lobe(_ln, 568.8, 46.9, 40.5) + 0.286 * _lobe(_ln, 530.9, 16.3, 31.1)
ZBAR = 1.217 * _lobe(_ln, 437.0, 11.8, 36.0) + 0.681 * _lobe(_ln, 459.0, 26.0, 13.8)
CMF = np.stack([XBAR, YBAR, ZBAR], -1)                       # [lam, 3]
XYZ2RGB = np.array([[3.2406, -1.5372, -0.4986], [-0.9689, 1.8758, 0.0415], [0.0557, -0.2040, 1.0570]])


def omega_geo_of_lambda(M_kg):
    """the geometric frequency w M of each wavelength on LAM, for a hole of mass M_kg"""
    unit = C ** 3 / (G * M_kg)                              # rad/s per unit geometric frequency
    return 2 * math.pi * C / LAM / unit


def xyz_of_radiance(I_of_w, M_kg):
    """XYZ of a radiance given per unit geometric frequency, I_of_w(w) with w an array, for a hole of mass M_kg"""
    w = omega_geo_of_lambda(M_kg)
    I_w = I_of_w(w)
    I_lam = I_w * w / LAM                                    # I_lam dlam = I_w dw, |dw/dlam| = w / lam
    dl = LAM[1] - LAM[0]
    return (I_lam[:, None] * CMF).sum(0) * dl


def rgb_of_xyz(xyz):
    return np.clip(XYZ2RGB @ xyz, 0, None)


def planck_w(w, T_geo):
    x = np.clip(w / T_geo, 1e-9, 700)
    return w ** 3 / (4 * math.pi ** 3 * np.expm1(x))


def blackbody_rgb(T_K, M_kg):
    """linear sRGB (with its absolute luminance in the common radiance units) of a blackbody at T_K"""
    T_geo = T_K / (MK_PER_TH / M_kg) * T_H_GEO              # T_K / T_H,K * T_H,geo
    return rgb_of_xyz(xyz_of_radiance(lambda w: planck_w(w, T_geo), M_kg))


class Hawking:
    def __init__(self, npz=None):
        npz = pathlib.Path(npz) if npz else pathlib.Path(__file__).resolve().parent / "greybody.npz"
        d = np.load(npz)
        self.w, self.S = d["w"], d["S"]
        self.w_peak = float(d["w_peak"])
        self.P, self.Ndot = float(d["P"]), float(d["Ndot"])
        self.d = d

    def greybody_ratio(self, w):
        """S(w) / (27 w^2): 1 in geometric optics, -> 0 at low frequency"""
        lo, hi = self.w[0], self.w[-1]
        ratio = self.S / (27 * self.w ** 2)
        out = np.interp(np.log(np.clip(w, lo, hi)), np.log(self.w), ratio)
        out = np.where(w < lo, ratio[0] * (w / lo) ** 2, out)      # Gamma_1 ~ w^4 below the grid: ratio ~ w^2
        out = np.where(w > hi, 1.0, out)
        return out

    def radiance_inf(self, w):
        return self.greybody_ratio(w) * planck_w(w, T_H_GEO)

    def rgb(self, M_kg, beta=1.0):
        """linear sRGB of the glow seen by a static observer with blueshift beta"""
        return rgb_of_xyz(xyz_of_radiance(lambda w: beta ** 3 * self.radiance_inf(w / beta), M_kg))

    def power_watts(self, M_kg):
        return self.P * HBAR * C ** 6 / (G ** 2 * M_kg ** 2)


def tonemap(rgb_lin):
    """luminance-preserving Reinhard, then the display gamma; rgb_lin [..., 3] linear, exposure already applied"""
    L = 0.2126 * rgb_lin[..., 0] + 0.7152 * rgb_lin[..., 1] + 0.0722 * rgb_lin[..., 2]
    Lt = L / (1.0 + L)
    scale = np.where(L > 1e-12, Lt / np.maximum(L, 1e-12), 0.0)
    out = np.clip(rgb_lin * scale[..., None], 0, 1)
    over = np.clip(L - 4.0, 0, None) / 8.0                      # the very bright bleach toward white, as eyes do
    out = out + (1.0 - out) * np.clip(over, 0, 1)[..., None]
    return np.clip(out, 0, 1) ** (1 / 2.2)


# ---- a 5x7 font ---------------------------------------------------------------------------------------------------
_F = {}
_F["A"] = ".###. #...# #...# ##### #...# #...# #...#"
_F["B"] = "####. #...# #...# ####. #...# #...# ####."
_F["C"] = ".###. #...# #.... #.... #.... #...# .###."
_F["D"] = "####. #...# #...# #...# #...# #...# ####."
_F["E"] = "##### #.... #.... ####. #.... #.... #####"
_F["F"] = "##### #.... #.... ####. #.... #.... #...."
_F["G"] = ".###. #...# #.... #.### #...# #...# .###."
_F["H"] = "#...# #...# #...# ##### #...# #...# #...#"
_F["I"] = ".###. ..#.. ..#.. ..#.. ..#.. ..#.. .###."
_F["J"] = "..### ...#. ...#. ...#. ...#. #..#. .##.."
_F["K"] = "#...# #..#. #.#.. ##... #.#.. #..#. #...#"
_F["L"] = "#.... #.... #.... #.... #.... #.... #####"
_F["M"] = "#...# ##.## #.#.# #.#.# #...# #...# #...#"
_F["N"] = "#...# ##..# #.#.# #..## #...# #...# #...#"
_F["O"] = ".###. #...# #...# #...# #...# #...# .###."
_F["P"] = "####. #...# #...# ####. #.... #.... #...."
_F["Q"] = ".###. #...# #...# #...# #.#.# #..#. .##.#"
_F["R"] = "####. #...# #...# ####. #.#.. #..#. #...#"
_F["S"] = ".#### #.... #.... .###. ....# ....# ####."
_F["T"] = "##### ..#.. ..#.. ..#.. ..#.. ..#.. ..#.."
_F["U"] = "#...# #...# #...# #...# #...# #...# .###."
_F["V"] = "#...# #...# #...# #...# #...# .#.#. ..#.."
_F["W"] = "#...# #...# #...# #.#.# #.#.# ##.## #...#"
_F["X"] = "#...# #...# .#.#. ..#.. .#.#. #...# #...#"
_F["Y"] = "#...# #...# .#.#. ..#.. ..#.. ..#.. ..#.."
_F["Z"] = "##### ....# ...#. ..#.. .#... #.... #####"
_F["a"] = "..... ..... .###. ....# .#### #...# .####"
_F["b"] = "#.... #.... ####. #...# #...# #...# ####."
_F["c"] = "..... ..... .###. #.... #.... #...# .###."
_F["d"] = "....# ....# .#### #...# #...# #...# .####"
_F["e"] = "..... ..... .###. #...# ##### #.... .###."
_F["f"] = "..##. .#..# .#... ###.. .#... .#... .#..."
_F["g"] = "..... .#### #...# #...# .#### ....# .###."
_F["h"] = "#.... #.... ####. #...# #...# #...# #...#"
_F["i"] = "..#.. ..... .##.. ..#.. ..#.. ..#.. .###."
_F["j"] = "...#. ..... ..##. ...#. ...#. #..#. .##.."
_F["k"] = "#.... #.... #..#. #.#.. ##... #.#.. #..#."
_F["l"] = ".##.. ..#.. ..#.. ..#.. ..#.. ..#.. .###."
_F["m"] = "..... ..... ##.#. #.#.# #.#.# #.#.# #...#"
_F["n"] = "..... ..... ####. #...# #...# #...# #...#"
_F["o"] = "..... ..... .###. #...# #...# #...# .###."
_F["p"] = "..... ####. #...# #...# ####. #.... #...."
_F["q"] = "..... .#### #...# #...# .#### ....# ....#"
_F["r"] = "..... ..... #.##. ##..# #.... #.... #...."
_F["s"] = "..... ..... .#### #.... .###. ....# ####."
_F["t"] = ".#... .#... ###.. .#... .#... .#..# ..##."
_F["u"] = "..... ..... #...# #...# #...# #..## .##.#"
_F["v"] = "..... ..... #...# #...# #...# .#.#. ..#.."
_F["w"] = "..... ..... #...# #...# #.#.# #.#.# .#.#."
_F["x"] = "..... ..... #...# .#.#. ..#.. .#.#. #...#"
_F["y"] = "..... #...# #...# #...# .#### ....# .###."
_F["z"] = "..... ..... ##### ...#. ..#.. .#... #####"
_F["0"] = ".###. #...# #..## #.#.# ##..# #...# .###."
_F["1"] = "..#.. .##.. ..#.. ..#.. ..#.. ..#.. .###."
_F["2"] = ".###. #...# ....# ...#. ..#.. .#... #####"
_F["3"] = "##### ...#. ..#.. ...#. ....# #...# .###."
_F["4"] = "...#. ..##. .#.#. #..#. ##### ...#. ...#."
_F["5"] = "##### #.... ####. ....# ....# #...# .###."
_F["6"] = "..##. .#... #.... ####. #...# #...# .###."
_F["7"] = "##### ....# ...#. ..#.. .#... .#... .#..."
_F["8"] = ".###. #...# #...# .###. #...# #...# .###."
_F["9"] = ".###. #...# #...# .#### ....# ...#. .##.."
_F[" "] = "..... ..... ..... ..... ..... ..... ....."
_F["."] = "..... ..... ..... ..... ..... .##.. .##.."
_F[","] = "..... ..... ..... ..... .##.. ..#.. .#..."
_F["="] = "..... ..... ##### ..... ##### ..... ....."
_F["-"] = "..... ..... ..... ##### ..... ..... ....."
_F["+"] = "..... ..#.. ..#.. ##### ..#.. ..#.. ....."
_F["/"] = "....# ....# ...#. ..#.. .#... #.... #...."
_F["("] = "..#.. .#... #.... #.... #.... .#... ..#.."
_F[")"] = "..#.. ...#. ....# ....# ....# ...#. ..#.."
_F[":"] = "..... .##.. .##.. ..... .##.. .##.. ....."
_F["^"] = "..#.. .#.#. #...# ..... ..... ..... ....."
_F["_"] = "..... ..... ..... ..... ..... ..... #####"
_F["%"] = "##... ##..# ...#. ..#.. .#... #..## ...##"
_F["|"] = "..#.. ..#.. ..#.. ..#.. ..#.. ..#.. ..#.."
_F["~"] = "..... ..... .#... #.#.# ...#. ..... ....."
_F["'"] = "..#.. ..#.. ..... ..... ..... ..... ....."
_F["<"] = "...#. ..#.. .#... #.... .#... ..#.. ...#."
_F[">"] = ".#... ..#.. ...#. ....# ...#. ..#.. .#..."
_F["*"] = "..... #.#.# .###. ##### .###. #.#.# ....."
_F["["] = ".###. .#... .#... .#... .#... .#... .###."
_F["]"] = ".###. ...#. ...#. ...#. ...#. ...#. .###."
GLYPH = {}
for ch, rows in _F.items():
    r = rows.split(" ")
    assert len(r) == 7 and all(len(x) == 5 for x in r), ch
    GLYPH[ch] = np.array([[c == "#" for c in row] for row in r])


def text_width(text, scale=2):
    return len(text) * 6 * scale


def draw_text(img, x, y, text, scale=2, color=(220, 220, 220)):
    """draw text into img (H, W, 3) uint8 at top-left (x, y); returns the x after the text"""
    H, W = img.shape[:2]
    col = np.array(color, np.uint8)
    for ch in text:
        g = GLYPH.get(ch, GLYPH["?"] if "?" in GLYPH else GLYPH[" "])
        big = np.kron(g, np.ones((scale, scale), bool))
        h, w = big.shape
        y0, x0 = y, x
        y1, x1 = min(H, y0 + h), min(W, x0 + w)
        if y1 > y0 and x1 > x0 and y0 >= 0 and x0 >= 0:
            sub = img[y0:y1, x0:x1]
            m = big[: y1 - y0, : x1 - x0]
            sub[m] = col
        x += 6 * scale
    return x


def write_ppm(path, img):
    path = pathlib.Path(path)
    with open(path, "wb") as f:
        f.write(f"P6 {img.shape[1]} {img.shape[0]} 255".encode() + bytes([10]))
        f.write(np.ascontiguousarray(img, dtype=np.uint8).tobytes())

"""The approach: a static (hovering) observer descends toward a Schwarzschild black hole that is radiating Hawking
photons into empty space (the Unruh state: outgoing modes thermal at T_H, nothing coming in but starlight), and
looks straight down and straight up (2026-09-08, the owner's question: can the never-seen Hawking glow be rendered
by going closer than anyone has looked).

    hawkapproach.py <outdir> [frames] [panel px] [T_H in kelvin] [r_start / M] [r_end / M]

Ray optics (honest for the short-wavelength tail of the spectrum, wrong at its peak, where the wavelength is 13 horizon
radii and the glow is a dipole blur; see hawkmode.py): every direction the observer looks is traced backwards, and it
either came from the horizon or from infinity. From the horizon it carries the Hawking glow: a uniform radiance,
thermal at T_H at infinity, blueshifted for the static observer at r by beta = 1 / sqrt(1 - 2M/r), the same in every
direction, so the glow is a uniform disc (or, inside r = 3M, a uniform sky with a hole in it) at the local temperature
T_H beta. From infinity it carries a faint star field, lensed, and blueshifted by the same beta. The border between
the two is the escape cone: from the outward direction, sky for angles psi < psi_e with sin(psi_e) = (3 sqrt3 M / r)
sqrt(1 - 2M/r) inside r = 3M, and psi_e = pi - asin(...) outside; every ray is also integrated numerically (the affine
null geodesic r'' = b^2 (1 - 3M/r) / r^3) for the lensing of the stars, and the two agree. The observer hovers: the
thrust needed is a = M / (r^2 sqrt(1 - 2M/r)), and as r -> 2M the glow's temperature T_H beta -> a / 2 pi, the Unruh
temperature of that thrust: at the horizon Hawking's radiation and Unruh's are one thing.

The glow's colour is the blackbody at T_H beta, integrated against the CIE colour-matching functions for a hole of
the mass that gives the chosen T_H (2000 K: 6.1e19 kg, a 91-nanometre horizon). Not the greybody-filtered spectrum
of greybody.npz: the filtering happens at the barrier at r = 3M, so a hoverer inside it sees the unfiltered thermal
flux, and the filter is the same wave effect that forbids this picture's sharp edge; the blackbody is the one spectrum
consistent with ray optics, and the real spectrum is in the chart. A camera's auto-exposure holds the glow at a fixed
display luminance (the caption counts the stops, which climb by 17 over the descent) so the colour is what changes:
orange-red at T_H, white near 5000 K, blue-white at 20,000 K. The stars are drawn at a fixed brightness with their
true colour shift (under that exposure they would fade out within a few M). Output: rgb24 frames
(<outdir>/frames.rgb), 1920 x 1040 (two 960 panels and a caption strip), and PPM stills at a few radii.
"""
import math
import pathlib
import sys
import time

import numpy as np

sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
import hawkcolor as hc

out = pathlib.Path(sys.argv[1]); out.mkdir(parents=True, exist_ok=True)
NF = int(sys.argv[2]) if len(sys.argv) > 2 else 240
P = int(sys.argv[3]) if len(sys.argv) > 3 else 960
TH_K = float(sys.argv[4]) if len(sys.argv) > 4 else 2000.0
R0 = float(sys.argv[5]) if len(sys.argv) > 5 else 40.0
R1 = float(sys.argv[6]) if len(sys.argv) > 6 else 2.02
FOV = math.radians(100.0)
HUD = 80
W, H = 2 * P, P + HUD
R_FAR = 800.0
SQ27 = 3 * math.sqrt(3)

t0 = time.time()
M_kg = hc.mass_for_TH(TH_K)
hk = hc.Hawking()
rs_m = hc.rs_metres(M_kg)
P_W = hk.power_watts(M_kg)
print(f"T_H = {TH_K:.0f} K -> M = {M_kg:.3e} kg, r_s = {rs_m * 1e9:.1f} nm, photon power {P_W:.3e} W; spectrum peak at w M = {hk.w_peak:.3f}")


def lum(rgb):
    return 0.2126 * rgb[..., 0] + 0.7152 * rgb[..., 1] + 0.0722 * rgb[..., 2]


glow1 = hc.blackbody_rgb(TH_K, M_kg)
TARGET_L = 0.8                       # the glow's linear luminance before the tone map (0.44 after it, 0.69 displayed)
EXPO0 = TARGET_L / lum(glow1)        # the exposure at the start; then a camera's auto-exposure holds the glow there
STAR_T = 3000.0 * (10000.0 / 3000.0) ** (np.arange(8) / 7.0)   # and the caption counts the stops
print(f"exposure: glow at T_H -> linear luminance {TARGET_L} (linear rgb of the glow at beta = 1: {np.round(glow1 / glow1.max(), 3)})")
# the stars are drawn at a fixed display brightness with their true colour shift: under the auto-exposure that holds
# the glow they would fade to nothing within a few M (the glow, in the Wien tail, brightens thousands of times faster)

# ---- the two panels' look directions (fixed) ------------------------------------------------------------------
tanx = math.tan(FOV / 2)
px = (np.arange(P) + 0.5) / P * 2 - 1
py = 1 - (np.arange(P) + 0.5) / P * 2
PX, PY = np.meshgrid(px, py)
panels = {}
for name, fz in (("down", -1.0), ("up", 1.0)):
    # the observer's frame: z is up (away from the hole), x right, y up-in-the-picture; looking down flips z and x
    nx = tanx * PX * (1.0 if fz > 0 else -1.0)
    ny = tanx * PY
    nz = np.full_like(nx, fz)
    nrm = np.sqrt(nx * nx + ny * ny + nz * nz)
    nx, ny, nz = nx / nrm, ny / nrm, nz / nrm
    psi = np.arccos(np.clip(nz, -1, 1))                              # angle from up
    s = np.maximum(np.sqrt(nx * nx + ny * ny), 1e-12)
    panels[name] = (psi.ravel(), (nx / s).ravel(), (ny / s).ravel())  # psi and the unit transverse direction e


def cone(r):
    B = (SQ27 / r) * math.sqrt(1 - 2 / r)
    return math.asin(min(B, 1.0)) if r < 3 else math.pi - math.asin(min(B, 1.0))


def trace(r_obs, psis):
    """null geodesics from a static observer at r_obs looking at angles psis from the outward radial direction.
    Returns (escaped mask, phi_infinity) with phi the position angle from the outward direction in the ray's plane."""
    n = psis.size
    b = r_obs * np.sin(psis) / math.sqrt(1 - 2 / r_obs)
    r = np.full(n, r_obs); pr = np.cos(psis); ph = np.zeros(n)
    state = np.zeros(n, np.int8)                                     # 0 flying, 1 escaped, 2 fell / lost
    phi_inf = np.zeros(n)
    active = np.arange(n)
    b2 = b * b
    steps = 0
    while active.size and steps < 60000:
        ra, pa, pha, bb = r[active], pr[active], ph[active], b2[active]
        h = np.clip(0.01 * ra, 0.004, 5.0)

        def f(rr, pp):
            return pp, bb * (1.0 / rr ** 3 - 3.0 / rr ** 4), np.sqrt(bb) / (rr * rr)

        k1r, k1p, k1f = f(ra, pa)
        k2r, k2p, k2f = f(ra + 0.5 * h * k1r, pa + 0.5 * h * k1p)
        k3r, k3p, k3f = f(ra + 0.5 * h * k2r, pa + 0.5 * h * k2p)
        k4r, k4p, k4f = f(ra + h * k3r, pa + h * k3p)
        rn = ra + h / 6 * (k1r + 2 * k2r + 2 * k3r + k4r)
        pn = pa + h / 6 * (k1p + 2 * k2p + 2 * k3p + k4p)
        phn = pha + h / 6 * (k1f + 2 * k2f + 2 * k3f + k4f)
        esc = rn > R_FAR
        fell = rn < 2.0005
        idx = active
        r[idx] = rn; pr[idx] = pn; ph[idx] = phn
        state[idx[esc]] = 1
        phi_inf[idx[esc]] = phn[esc] + np.arcsin(np.clip(np.sqrt(bb[esc]) / rn[esc], -1, 1))
        state[idx[fell & ~esc]] = 2
        active = idx[~(esc | fell)]
        steps += 1
    state[active] = 2
    return state == 1, phi_inf


def deflection_table(r_obs, psi_e):
    qa = np.linspace(0, 1, 1400, endpoint=False) * psi_e
    qb = psi_e - np.exp(np.linspace(math.log(psi_e / 1400), math.log(1e-8), 1400))
    psis = np.unique(np.concatenate([qa, qb]))
    psis = psis[(psis > 1e-5) & (psis < psi_e)]
    escm, phi_inf = trace(r_obs, psis)
    # the nearly radial rays are straight
    psis = np.concatenate([[0.0], psis[escm]]); phi_inf = np.concatenate([[0.0], phi_inf[escm]])
    return psis, phi_inf, int((~escm).sum())


def stars(ex, ey, ez):
    v0 = np.floor((ex + 1.0) * 400.0).astype(np.int64)
    v1 = np.floor((ey + 1.0) * 400.0).astype(np.int64)
    v2 = np.floor((ez + 1.0) * 400.0).astype(np.int64)
    h = (v0 * 73856093 ^ v1 * 19349663 ^ v2 * 83492791) & 0xFFFFFF
    h2 = (v0 * 83492791 ^ v1 * 73856093 ^ v2 * 19349663) & 0xFFFF
    on = h < int(0xFFFFFF * 0.004)
    bright = 0.25 + 0.75 * ((h2 >> 3) & 0xFF) / 255.0
    cls = h2 & 7
    return on, bright, cls


# ---- the frames ------------------------------------------------------------------------------------------------
rr = 2.0 + (R0 - 2.0) * ((R1 - 2.0) / (R0 - 2.0)) ** (np.arange(NF) / max(NF - 1, 1))
stills = {int(np.argmin(np.abs(rr - v))): v for v in (R0, 10.0, 4.0, 2.5, 2.1, R1)}
fh = open(out / "frames.rgb", "wb")
checked = False
for n in range(NF):
    r = float(rr[n])
    beta = 1.0 / math.sqrt(1 - 2 / r)
    T_loc = TH_K * beta
    psi_e = cone(r)
    accel = 1.0 / (r * r * math.sqrt(1 - 2 / r))                   # in units of c^4 / (G M)
    accel_g = accel * hc.C ** 4 / (hc.G * M_kg) / 9.80665
    psis, phis, lost = deflection_table(r, psi_e)
    if not checked:
        # the analytic cone against the tracer, once: a ray just inside the cone escapes, just outside falls
        e_in, _ = trace(r, np.array([psi_e - 0.02]))
        e_out, _ = trace(r, np.array([min(psi_e + 0.02, math.pi - 1e-4)]))
        print(f"cone check at r = {r:.2f}: psi_e = {math.degrees(psi_e):.2f} deg; inside escapes: {bool(e_in[0])}, outside escapes: {bool(e_out[0])}")
        checked = True
    glow_lin = hc.blackbody_rgb(T_loc, M_kg)
    EXPO = TARGET_L / lum(glow_lin)
    stops = math.log2(EXPO0 / EXPO)
    glow = glow_lin * EXPO
    star_rgb = []
    for T in STAR_T:
        c = hc.blackbody_rgb(float(T) * beta, M_kg)
        star_rgb.append(c / max(lum(c), 1e-30) * 0.9)
    star_rgb = np.stack(star_rgb)
    frame = np.zeros((H, W, 3), np.uint8)
    for k, name in enumerate(("down", "up")):
        psi, ex, ey = panels[name]
        sky = psi < psi_e
        rgb = np.tile(glow[None, :], (psi.size, 1))
        if sky.any():
            phi_inf = np.interp(psi[sky], psis, phis)
            c, s = np.cos(phi_inf), np.sin(phi_inf)
            on, bright, cls = stars(s * ex[sky], s * ey[sky], c)
            col = star_rgb[cls] * (bright * on)[:, None]
            rgb[sky] = col
        img = hc.tonemap(rgb.reshape(P, P, 3))
        frame[:P, k * P:(k + 1) * P] = (img * 255 + 0.5).astype(np.uint8)
    frame[P - 1:P + 1, :] = 40
    frame[:P, P - 1:P + 1] = 40
    hc.draw_text(frame, 16, 14, "looking down at the hole", 2, (200, 200, 200))
    hc.draw_text(frame, P + 16, 14, "looking up, away from it", 2, (200, 200, 200))
    line1 = f"hovering at r = {r:6.3f} M    glow T = {beta:5.2f} T_H = {T_loc:6.0f} K    sky cone {math.degrees(psi_e):5.1f} deg    exposure {-stops:+5.1f} stops"
    line2 = f"thrust to hover a = {accel:8.4f} c^4/GM = {accel_g:8.2e} g    T_H = {TH_K:.0f} K, M = {M_kg:.2e} kg, horizon {rs_m * 1e9:.0f} nm, {P_W * 1e9:.0f} nW    stars: fixed brightness, true colour shift"
    hc.draw_text(frame, 16, P + 14, line1, 2, (230, 230, 230))
    hc.draw_text(frame, 16, P + 46, line2, 2, (150, 150, 150))
    fh.write(frame.tobytes())
    if n in stills:
        hc.write_ppm(out / f"still_r{stills[n]:.2f}.ppm", frame)
    if n % 24 == 0 or n == NF - 1:
        print(f"  frame {n:3d}: r = {r:6.3f}, beta = {beta:5.2f}, T = {T_loc:6.0f} K, cone {math.degrees(psi_e):5.1f} deg, table {psis.size} rays ({lost} lost) ({time.time() - t0:.0f} s)", flush=True)
fh.close()
print(f"{NF} frames of {W}x{H} -> {out / 'frames.rgb'} ({time.time() - t0:.0f} s)")

"""The Hawking mode itself, at the wavelength it is actually emitted at (2026-09-08).

    hawkmode.py <outdir> [frames] [panel width] [panel height] [pixels per M]

Ray optics cannot picture Hawking radiation at its peak: the photon spectrum peaks at w M = 0.243 (greybody.npz), a
wavelength of 25.8 M = 12.9 horizon radii, so the emitter is far smaller than its own light and only the l = 1 partial
wave gets out (98 % of the power). What CAN be drawn exactly is the mode function: the electromagnetic "up" mode of
the Regge-Wheeler equation at that frequency, born at the horizon with unit outgoing amplitude, partly transmitted
through the potential barrier that peaks at r = 3M (Gamma_1 = 0.415), partly reflected back in. Left panel l = m = 1,
right panel l = m = 2 (Gamma_2 = 0.0004: trapped), both in the equatorial plane, the field Re[psi(r*) exp(i (m phi -
w t))] as a signed colour (warm positive, cool negative), the flux-normalised amplitude psi rather than the 1/r field
so the outgoing wave keeps its brightness. Near the horizon the wave's crests pile up (r* -> -infinity: the same
outgoing wave, infinitely compressed in r, the trans-Planckian side of Hawking's derivation) and peel off it at the
coordinate speed 1 - 2M/r. The horizon is the black disc, the barrier's peak the dashed ring. Output rgb24 frames.
"""
import math
import pathlib
import sys
import time

import numpy as np

sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
import hawkcolor as hc

out = pathlib.Path(sys.argv[1]); out.mkdir(parents=True, exist_ok=True)
NF = int(sys.argv[2]) if len(sys.argv) > 2 else 240
PW = int(sys.argv[3]) if len(sys.argv) > 3 else 960
PH = int(sys.argv[4]) if len(sys.argv) > 4 else 1080
PPM = float(sys.argv[5]) if len(sys.argv) > 5 else 16.0
PERIODS = 3.0

t0 = time.time()
d = np.load(pathlib.Path(__file__).resolve().parent / "greybody.npz")
w = float(d["w_peak"])
modes = {1: (d["rs1"], d["psi1"], float(d["gamma1"])), 2: (d["rs2"], d["psi2"], float(d["gamma2"]))}
lam = 2 * math.pi / w
print(f"w = {w:.4f} / M, wavelength {lam:.1f} M = {lam / 2:.1f} r_s; Gamma_1 = {modes[1][2]:.4f}, Gamma_2 = {modes[2][2]:.5f}")

xs = (np.arange(PW) + 0.5 - PW / 2) / PPM
ys = -(np.arange(PH) + 0.5 - PH / 2) / PPM
X, Y = np.meshgrid(xs, ys)
R = np.hypot(X, Y)
PHI = np.arctan2(Y, X)
outside = R > 2.0
xh = np.maximum(R - 2.0, 1e-12)
RS = 2.0 + xh + 2.0 * np.log(xh / 2.0)
fields = {}
for l, (rs, psi, g) in modes.items():
    A = np.interp(RS, rs, psi.real)
    B = np.interp(RS, rs, psi.imag)
    A[~outside] = 0; B[~outside] = 0
    fields[l] = (A, B)
vmax = max(np.sqrt(A * A + B * B).max() for A, B in fields.values())
print(f"amplitude range: max |psi| = {vmax:.3f}; l = 1 outside the barrier |psi| = {np.sqrt(fields[1][0] ** 2 + fields[1][1] ** 2)[(R > 6) & (R < 8)].mean():.3f}, l = 2 there {np.sqrt(fields[2][0] ** 2 + fields[2][1] ** 2)[(R > 6) & (R < 8)].mean():.4f}")

WARM = np.array([1.0, 0.62, 0.22])
COOL = np.array([0.25, 0.55, 1.0])
horizon_ring = np.abs(R - 2.0) < 0.75 / PPM
barrier = (np.abs(R - 3.0) < 0.6 / PPM) & ((np.floor(PHI / (math.pi / 18)).astype(int) % 2) == 0)
period = 2 * math.pi / w
fh = open(out / "frames.rgb", "wb")
for n in range(NF):
    t = n / NF * PERIODS * period
    frame = np.zeros((PH, 2 * PW, 3), np.uint8)
    for k, l in enumerate((1, 2)):
        A, B = fields[l]
        ph = l * PHI - w * t
        v = (A * np.cos(ph) - B * np.sin(ph)) / vmax
        a = np.abs(v) ** 0.75
        rgb = np.where((v > 0)[..., None], WARM, COOL) * a[..., None]
        rgb[~outside] = 0
        rgb[horizon_ring] = 0.85
        rgb[barrier] = np.maximum(rgb[barrier], 0.45)
        img = (np.clip(rgb, 0, 1) ** (1 / 2.2) * 255 + 0.5).astype(np.uint8)
        frame[:, k * PW:(k + 1) * PW] = img
    frame[:, PW - 1:PW + 1] = 40
    # the wavelength bar and captions
    x0, y0 = 24, 110
    frame[y0 - 1:y0 + 2, x0:x0 + int(lam * PPM)] = 220
    frame[y0 - 8:y0 + 9, x0:x0 + 2] = 220; frame[y0 - 8:y0 + 9, x0 + int(lam * PPM) - 2:x0 + int(lam * PPM)] = 220
    hc.draw_text(frame, x0, y0 - 34, f"one wavelength: {lam:.1f} M = {lam / 2:.1f} horizon radii", 2, (220, 220, 220))
    hc.draw_text(frame, 24, 14, "the Hawking photon mode at its spectrum's peak, w = 0.243 / M", 2, (200, 200, 200))
    hc.draw_text(frame, 24, 40, "flux-normalised amplitude psi, equatorial plane, time e^(-i w t)", 2, (150, 150, 150))
    hc.draw_text(frame, PW + 24, 14, "the next multipole at the same frequency", 2, (200, 200, 200))
    hc.draw_text(frame, 24, PH - 60, f"l = 1, m = 1    Gamma = {modes[1][2]:.3f}: the dipole gets out (98% of the power)", 2, (230, 230, 230))
    hc.draw_text(frame, 24, PH - 32, "black disc: horizon r = 2 M    dashed: barrier peak r = 3 M", 2, (150, 150, 150))
    hc.draw_text(frame, PW + 24, PH - 60, f"l = 2, m = 2    Gamma = {modes[2][2]:.4f}: the quadrupole is trapped", 2, (230, 230, 230))
    hc.draw_text(frame, PW + 24, PH - 32, "a near-standing wave inside the barrier; a few percent leak out", 2, (150, 150, 150))
    fh.write(frame.tobytes())
    if n == NF // 2:
        hc.write_ppm(out / "still_mid.ppm", frame)
    if n % 60 == 0:
        print(f"  frame {n} ({time.time() - t0:.0f} s)", flush=True)
fh.close()
print(f"{NF} frames of {2 * PW}x{PH} -> {out / 'frames.rgb'} ({time.time() - t0:.0f} s)")

"""The chart: photon greybody factors and the Hawking photon spectrum against the blackbody it is usually drawn as,
plus the mass ladder printed for the record (2026-09-08). numpy only, drawn by hand; writes <out.ppm>.

    hawkchart.py <out.ppm>
"""
import math
import pathlib
import sys

import numpy as np

sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
import hawkcolor as hc

outp = pathlib.Path(sys.argv[1])
d = np.load(pathlib.Path(__file__).resolve().parent / "greybody.npz")
w, Gamma, dEdw, dEdw_geo = d["w"], d["Gamma"], d["dEdw"], d["dEdw_geo"]
P, Ndot, w_pk = float(d["P"]), float(d["Ndot"]), float(d["w_peak"])
T_H = 1 / (8 * math.pi)

W, H = 1600, 1100
img = np.zeros((H, W, 3), np.uint8)
img[:] = 12


class Axes:
    def __init__(self, box, xr, yr, logx=False, logy=False):
        self.x0, self.y0, self.x1, self.y1 = box
        self.xr, self.yr, self.logx, self.logy = xr, yr, logx, logy

    def _tx(self, x):
        a, b = self.xr
        if self.logx:
            x, a, b = np.log(x), math.log(a), math.log(b)
        return self.x0 + (x - a) / (b - a) * (self.x1 - self.x0)

    def _ty(self, y):
        a, b = self.yr
        if self.logy:
            y, a, b = np.log(np.maximum(y, 1e-300)), math.log(a), math.log(b)
        return self.y1 - (y - a) / (b - a) * (self.y1 - self.y0)

    def frame(self):
        img[self.y0:self.y1 + 1, self.x0:self.x0 + 2] = 120
        img[self.y0:self.y1 + 1, self.x1 - 1:self.x1 + 1] = 120
        img[self.y0:self.y0 + 2, self.x0:self.x1 + 1] = 120
        img[self.y1 - 1:self.y1 + 1, self.x0:self.x1 + 1] = 120

    def line(self, xs, ys, color, width=2):
        px, py = self._tx(np.asarray(xs, float)), self._ty(np.asarray(ys, float))
        ok = np.isfinite(px) & np.isfinite(py)
        px, py = px[ok], py[ok]
        for i in range(len(px) - 1):
            n = int(max(2, math.hypot(px[i + 1] - px[i], py[i + 1] - py[i]) * 1.5))
            xx = np.linspace(px[i], px[i + 1], n); yy = np.linspace(py[i], py[i + 1], n)
            for dx in range(width):
                for dy in range(width):
                    xi = (xx + dx - width // 2).astype(int); yi = (yy + dy - width // 2).astype(int)
                    m = (xi >= self.x0) & (xi <= self.x1) & (yi >= self.y0) & (yi <= self.y1)
                    img[yi[m], xi[m]] = color

    def vband(self, xa, xb, color):
        a, b = int(self._tx(xa)), int(self._tx(xb))
        img[self.y0 + 2:self.y1 - 1, a:b] = color

    def xtick(self, x, label):
        px = int(self._tx(x))
        img[self.y1 - 8:self.y1, px:px + 2] = 160
        hc.draw_text(img, px - hc.text_width(label, 2) // 2, self.y1 + 8, label, 2, (190, 190, 190))

    def ytick(self, y, label):
        py = int(self._ty(y))
        img[py:py + 2, self.x0:self.x0 + 8] = 160
        hc.draw_text(img, self.x0 - hc.text_width(label, 2) - 8, py - 7, label, 2, (190, 190, 190))

    def text(self, x, y, s, color=(220, 220, 220), scale=2):
        hc.draw_text(img, int(self._tx(x)), int(self._ty(y)), s, scale, color)


hc.draw_text(img, 40, 24, "Hawking radiation of a Schwarzschild black hole in photons  (computed this session; Page 1976 reproduced)", 2, (240, 240, 240))

# top: greybody factors
ax = Axes((160, 90, 1540, 470), (0.01, 1.5), (1e-8, 1.5), logx=True, logy=True)
ax.frame()
cols = [(255, 170, 60), (90, 200, 255), (150, 255, 150), (255, 120, 200), (200, 200, 120)]
for l in range(1, 6):
    ax.line(w, Gamma[l], cols[l - 1], 3)
    i = int(np.argmin(np.abs(np.log(Gamma[l]) - math.log(1e-3))))
    ax.text(w[i] * 1.12, 1e-3, f"l = {l}", cols[l - 1])
for v in (0.01, 0.03, 0.1, 0.3, 1.0):
    ax.xtick(v, f"{v:g}")
for v in (1e-8, 1e-6, 1e-4, 1e-2, 1.0):
    ax.ytick(v, f"{v:.0e}" if v < 0.5 else "1")
hc.draw_text(img, 170, 100, "greybody factor Gamma_l(w): the fraction of a horizon-born wave that reaches infinity", 2, (220, 220, 220))
hc.draw_text(img, 700, 480 + 30, "frequency w in units of 1/M   (T_H = 1/8 pi M = 0.040/M)", 2, (190, 190, 190))

# bottom: the spectrum
ymax = float(dEdw_geo.max()) * 1.08
ax2 = Axes((160, 600, 1540, 1000), (0.0, 1.2), (0.0, ymax))
ax2.frame()
M_film = hc.mass_for_TH(2000.0)
unit = hc.C ** 3 / (hc.G * M_film)
wa, wb = 2 * math.pi * hc.C / 780e-9 / unit, 2 * math.pi * hc.C / 380e-9 / unit
ax2.vband(wa, wb, (40, 40, 52))
ax2.line(w, dEdw_geo, (110, 110, 110), 3)
ax2.line(w, dEdw, (255, 170, 60), 4)
ipk = int(np.argmax(dEdw))
ax2.line([w_pk, w_pk], [0, dEdw[ipk]], (255, 170, 60), 1)
for v in (0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2):
    ax2.xtick(v, f"{v:g}")
ax2.ytick(0.0, "0"); ax2.ytick(ymax / 1.08, f"{ymax / 1.08:.1e}")
hc.draw_text(img, 170, 610, "power spectrum dE/dt dw, in units of hbar c^6 / G^2 M^2 per unit w M", 2, (220, 220, 220))
hc.draw_text(img, 700, 1010 + 30, "frequency w in units of 1/M", 2, (190, 190, 190))
ax2.text(0.42, ymax * 0.92, "grey: blackbody at T_H over the capture area 27 pi M^2  (P = 1.40e-4)", (170, 170, 170))
ax2.text(0.42, ymax * 0.82, f"orange: with the greybody factors  (P = {P:.3e}; Page 1976: 3.36e-5)", (255, 170, 60))
ax2.text(0.42, ymax * 0.72, f"peak w = {w_pk:.3f}/M = {w_pk / T_H:.1f} T_H (blackbody 2.8 T_H); wavelength {2 * math.pi / w_pk:.1f} M = {math.pi / w_pk:.1f} r_s", (255, 170, 60))
ax2.text(0.42, ymax * 0.62, "98% of the power in l = 1: at its peak the hole is a pure dipole", (255, 170, 60))
ax2.text(wa + 0.005, ymax * 0.30, "visible band (380-780 nm)", (140, 140, 170))
ax2.text(wa + 0.005, ymax * 0.24, "for the film's hole: T_H = 2000 K", (140, 140, 170))
ax2.text(wa + 0.005, ymax * 0.18, "M = 6.1e19 kg, horizon 91 nm", (140, 140, 170))
hc.write_ppm(outp, img)
print(f"wrote {outp}")

# the mass ladder, for the record
print()
print("mass ladder (Schwarzschild, photons only; P from this session's 3.364e-5 hbar c^6 / G^2 M^2; peak wavelength 25.8 GM/c^2):")
print(f"{'M (kg)':>10} {'T_H (K)':>10} {'horizon':>10} {'peak lambda':>12} {'P (W)':>10} {'at 1 m (W/m^2)':>15}  note")
M_sun, M_earth, M_moon = 1.989e30, 5.972e24, 7.35e22
rows = [(1e17, ""), (1e18, ""), (1e19, ""), (hc.mass_for_TH(6000.0), "T_H = 6000 K: the Sun's colour"),
        (hc.mass_for_TH(2000.0), "the film's hole"), (1e21, ""), (1e22, ""), (hc.mass_for_TH(2.725), "T_H = the CMB: heavier holes absorb more than they emit"),
        (M_moon, "the Moon"), (M_earth, "the Earth"), (M_sun, "the Sun"), (4e6 * M_sun, "Sgr A*")]
hk = hc.Hawking()
for M, note in rows:
    T = hc.MK_PER_TH / M
    rs = hc.rs_metres(M)
    lam = 25.8 * hc.G * M / hc.C ** 2
    Pw = hk.power_watts(M)
    def fm(x, u=("m", "mm", "um", "nm", "pm", "fm", "am", "zm")):
        if x >= 1e3:
            return f"{x / 1e3:.3g} km"
        for i, s in enumerate(u):
            if x >= 10 ** (-3 * i) or i == len(u) - 1:
                return f"{x / 10 ** (-3 * i):.3g} {s}"
    print(f"{M:10.2e} {T:10.3g} {fm(rs):>10} {fm(lam):>12} {Pw:10.2e} {Pw / (4 * math.pi):15.2e}  {note}")

#!/bin/bash
# Hawking radiation, rendered as far as the mathematics honestly allows (2026-09-08): the hovering approach film,
# the mode film, the chart; encode and copy to hot-drops. CPU only; run alone.
set -u
FFDIR="${FFDIR:-/path/to/ffmpeg}"          # the directory holding ffmpeg.exe / ffprobe.exe
export PATH="$FFDIR:$PATH"
G=$(pwd); HD=$G/out; mkdir -p "$HD"
cd "$G"
echo "########## HAWKRENDER start $(date +%T)"
python hawkapproach.py hawk_approach 240 960 2000 40 2.02 2>&1 | tee hawk_approach.log
python hawkmode.py hawk_mode 240 960 1080 16 2>&1 | tee hawk_mode.log
python hawkchart.py hawk_chart.ppm 2>&1 | tee hawk_chart.log
echo "########## encode $(date +%T)"
ffmpeg.exe -y -hide_banner -loglevel error -f rawvideo -pix_fmt rgb24 -s 1920x1040 -r 24 -i hawk_approach/frames.rgb -vf format=yuv420p -c:v h264_mf -b:v 20M -movflags +faststart "$HD/hawking-approach-hover-40M-to-2.02M-10s.mp4" && echo "  approach: $(stat -c %s "$HD/hawking-approach-hover-40M-to-2.02M-10s.mp4") bytes"
ffmpeg.exe -y -hide_banner -loglevel error -f rawvideo -pix_fmt rgb24 -s 1920x1080 -r 24 -i hawk_mode/frames.rgb -vf format=yuv420p -c:v h264_mf -b:v 20M -movflags +faststart "$HD/hawking-mode-dipole-quadrupole-10s.mp4" && echo "  mode: $(stat -c %s "$HD/hawking-mode-dipole-quadrupole-10s.mp4") bytes"
for f in hawk_approach/still_r*.ppm; do
  b=$(basename "$f" .ppm); r=${b#still_r}
  ffmpeg.exe -y -hide_banner -loglevel error -i "$f" "$HD/hawking-approach-r${r}M.png" && echo "  still r = $r"
done
ffmpeg.exe -y -hide_banner -loglevel error -i hawk_mode/still_mid.ppm "$HD/hawking-mode-frame120.png" && echo "  mode still"
ffmpeg.exe -y -hide_banner -loglevel error -i hawk_chart.ppm "$HD/hawking-spectrum-chart.png" && echo "  chart"
ls -la "$HD" | grep hawking
echo "HAWKRENDER DONE $(date +%T)"

l = 1: Gamma at w = 0.010: 7.622e-08, at w = 0.124: 0.0064, at w = 1.50: 1.000150; flux conservation max |1 - (|A|^2 - |B|^2)| = 1.5e-04 (3 s)
l = 2: Gamma at w = 0.010: 7.585e-14, at w = 0.124: 0.0000, at w = 1.50: 1.000150; flux conservation max |1 - (|A|^2 - |B|^2)| = 5.9e-03 (5 s)
l = 3: Gamma at w = 0.010: 4.434e-20, at w = 0.124: 0.0000, at w = 1.50: 1.000150; flux conservation max |1 - (|A|^2 - |B|^2)| = 8.2e+03 (8 s)
l = 4: Gamma at w = 0.010: 1.554e-26, at w = 0.124: 0.0000, at w = 1.50: 1.000150; flux conservation max |1 - (|A|^2 - |B|^2)| = 2.1e+10 (11 s)
l = 5: Gamma at w = 0.010: 3.688e-33, at w = 0.124: 0.0000, at w = 1.50: 1.000149; flux conservation max |1 - (|A|^2 - |B|^2)| = 1.8e+17 (14 s)
l = 6: Gamma at w = 0.010: 6.139e-40, at w = 0.124: 0.0000, at w = 1.50: 0.999874; flux conservation max |1 - (|A|^2 - |B|^2)| = 1.5e+24 (16 s)
l = 7: Gamma at w = 0.010: 7.414e-47, at w = 0.124: 0.0000, at w = 1.50: 0.887382; flux conservation max |1 - (|A|^2 - |B|^2)| = 6.3e+30 (19 s)
l = 8: Gamma at w = 0.010: 7.492e-54, at w = 0.124: 0.0000, at w = 1.50: 0.013485; flux conservation max |1 - (|A|^2 - |B|^2)| = 8.5e+37 (22 s)
low-w slope of Gamma_1: 4.143 (theory 4)
  sum (2l+1) Gamma at w = 0.30: 2.511  vs geometric 27 w^2 = 2.448  (ratio 1.026)
  sum (2l+1) Gamma at w = 0.60: 8.928  vs geometric 27 w^2 = 9.858  (ratio 0.906)
  sum (2l+1) Gamma at w = 1.00: 25.695  vs geometric 27 w^2 = 26.957  (ratio 0.953)
  sum (2l+1) Gamma at w = 1.50: 61.544  vs geometric 27 w^2 = 60.750  (ratio 1.013)
photon power P M^2 = 3.364e-05  (Page 1976: 3.36e-5; geometric-optics blackbody 1.398e-04, exact 1.399e-04; ratio P/P_geo = 0.240)
photon rate  N M   = 1.480e-04 per unit time M;  mean photon energy 0.2273 = 5.71 T_H
energy spectrum peak at w M = 0.243 = 6.12 T_H  (blackbody peak 2.82 T_H = 0.112);  wavelength 2 pi / w = 25.8 M = 12.9 r_s
number spectrum peak at w M = 0.230
share of the power by multipole: l=1: 98.0%, l=2: 2.0%, l=3: 0.0%, l=4: 0.0%
mode l = 1 at w = 0.243: Gamma = 0.4148, |R|^2 = 0.5852, 2433 points kept, r* in [-31.6, 90.0] (24 s)
mode l = 2 at w = 0.243: Gamma = 0.0004, |R|^2 = 0.9996, 2433 points kept, r* in [-31.6, 90.0] (26 s)
wrote greybody.npz (26 s)

jp a=0.9 eps3=3.0 horizon=1.4359 isco=1.4646 hits=3498121 fell=150703 escaped=4645576 g_max=1.320
traced 8294400 rays in 6238 s: 3498121 hit the disc, 150703 fell in, 4645576 escaped (3840x2160, 240 frames, 4 bands)