M87 Black Hole

I asked Chatgpt 5 to write a Python script that replicates the First Image of a Black Hole, but using all 24-bit rgb colors. The result is not at all what I expected, but I like it a lot and so I decided to publish it.

It even gave a little explanation as to why this image is accurate. I like how it mentions the number 42. Who knows, maybe it stumbled upon something meaningful.

#!/usr/bin/env python3
"""
M87* all-RGB image (4096×4096), scientifically-informed crescent.

What this does
--------------
1) Builds a physically-motivated intensity map I(x,y):
   - Thin Gaussian ring ("photon ring" scale) with a strong SOUTHERN brightness
     asymmetry (relativistic beaming) and a deep central brightness depression.
   - Parameters guided by EHT measurements (ring diameter ~42 μas, thin ring).
     Refs: EHT 2019 reports; PRIMO 2023 sharper ring; 2024 multi-epoch ring
     with southerly brightness. See citations in the chat message.

2) Uses ALL 24-bit RGB colors exactly once:
   - Total pixels = 4096×4096 = 16,777,216.
   - Selects exactly as many “ring pixels” as needed by a geometric ring mask.
   - Assigns the most chromatic colors to the ring (ordered by hue along angle).
   - Assigns the remaining colors to all other pixels, dark-to-bright to match I.

Design choices
--------------
* Ring geometry:
  - r0_frac ~ 0.29 (radius as fraction of image size), thin width.
  - Southern bright sector via a Doppler-like factor A(θ) = (1+k cos(θ-θ0))^p,
    with θ0 = south (−π/2) in image coordinates (y down).
* Exact-use guarantee:
  - Colors are split into two non-overlapping pools (ring vs non-ring).
  - Each pool is strictly permutation-assigned to its pixel set.
  - A final uniqueness assertion checks we used every RGB exactly once.

Runtime & memory
----------------
This is heavy by design (tens to hundreds of MB of RAM and big sorts).
For lower RAM, you can chunk some steps, but this single-pass version
is easiest to verify for exact color accounting.

Outputs
-------
- PNG named: M87_allrgb_.png
- Console prints: counts & verification.
"""

import numpy as np
from PIL import Image
import time
import math

# ------------------------- parameters you might tweak -------------------------

SIZE = 4096  # image width=height; must be 4096 to match 16,777,216 colors
SAVE_TIMESTAMP = time.strftime("%Y%m%d_%H%M%S")
OUTFILE = f"M87_allrgb_{SAVE_TIMESTAMP}.png"

# Ring geometry (fraction of SIZE)
r0_frac    = 0.29   # ring radius
ring_band  = 96.0   # half-thickness of the ring band in pixels (visual ring thickness)
sigma_px   = 28.0   # Gaussian radial sigma controlling the ring profile thickness

# Beaming (brightness asymmetry, brightest toward south)
theta0     = -math.pi/2       # south (in image coords: +y down)
k_beam     = 0.65             # 0..1 strength of cosine beaming
p_beam     = 2.2              # raises contrast of beaming

# Background/contrast shaping
bg_envelope_strength = 0.25   # faint background glow
shadow_sharpness     = 8.0    # higher -> steeper dark center

# -----------------------------------------------------------------------------

# --- 1) Build coordinate grid and physically-motivated intensity I(x,y) ---

H = W = SIZE
yy, xx = np.indices((H, W), dtype=np.float32)
cx = (W - 1) * 0.5
cy = (H - 1) * 0.5

# Pixel-space radius & angle (θ measured from +x to -y for "image up")
dx = xx - cx
dy = yy - cy
r = np.hypot(dx, dy)
theta = np.arctan2(-dy, dx)  # so that -pi/2 points downward (south) in the image

r0 = r0_frac * SIZE

# Radial "thin ring": Gaussian around r0
I_rad = np.exp(-0.5 * ((r - r0) / sigma_px) ** 2)

# Doppler-like beaming: brightest near theta0 (south)
beam = 1.0 + k_beam * np.cos(theta - theta0)
beam = np.clip(beam, 0.0, None) ** p_beam

# Deep central depression (shadow) via a sharp logistic around r ~ r0 - a few sigmas
# This suppresses the interior brightness by >10:1 relative to ring.
inner_edge = r0 - 2.2 * sigma_px
shadow = 1.0 / (1.0 + np.exp(-(r - inner_edge) / (sigma_px / shadow_sharpness)))

# Faint background envelope to avoid a totally flat outer region
bg = bg_envelope_strength * np.exp(- (r / (0.42 * SIZE)) ** 2)

# Combine and normalize
I = (I_rad * beam * shadow) + bg
I -= I.min()
I /= (I.max() + 1e-12)

# --- ring mask: geometric band around r0 (uses band, not intensity threshold) ---
ring_mask = (np.abs(r - r0) <= ring_band) & (I > 0.02)
ring_idx = np.flatnonzero(ring_mask)
non_idx  = np.flatnonzero(~ring_mask)

N_total = SIZE * SIZE
N_ring  = ring_idx.size
N_non   = non_idx.size
assert N_ring + N_non == N_total

# Order ring pixels around azimuth, starting at theta0 (south) and going CCW
phi = (theta - theta0) % (2.0 * np.pi)
phi_ring = phi.flat[ring_idx]
order_ring_pixels = ring_idx[np.argsort(phi_ring, kind='mergesort')]

# Order non-ring pixels by intensity (dark -> bright)
I_non = I.flat[non_idx]
order_non_pixels = non_idx[np.argsort(I_non, kind='mergesort')]

# --- 2) Generate ALL 24-bit colors (no repeats), then split into two pools ---

# Colors as uint8 channels via bit tricks; stable and memory-light(ish)
N = 1 << 24  # 16,777,216
ids = np.arange(N, dtype=np.uint32)
R = ((ids >> 16) & 0xFF).astype(np.uint8)
G = ((ids >>  8) & 0xFF).astype(np.uint8)
B = ( ids        & 0xFF).astype(np.uint8)

# "Chromaticity" metric to pick the most colorful colors for the ring:
# Using max(R,G,B) - min(R,G,B) (channel spread). This equals (S*V)*255 in sRGB scale.
maxc = np.maximum(np.maximum(R, G), B).astype(np.int16)
minc = np.minimum(np.minimum(R, G), B).astype(np.int16)
spread = (maxc - minc).astype(np.int16)   # 0..255

# Select EXACTLY N_ring "most chromatic" colors for the ring
# (Tie-breaking handled by argpartition+argsort for determinism.)
if N_ring > 0:
    part_idx = np.argpartition(spread, -N_ring)[-N_ring:]
    # Now order those ring colors by hue for a smooth hue cycle around the ring
    # Compute hue ONLY for the selected subset to save memory/time.
    r_f = R[part_idx].astype(np.float32) / 255.0
    g_f = G[part_idx].astype(np.float32) / 255.0
    b_f = B[part_idx].astype(np.float32) / 255.0

    mx = np.maximum(np.maximum(r_f, g_f), b_f)
    mn = np.minimum(np.minimum(r_f, g_f), b_f)
    df = mx - mn
    # Hue in [0,1)
    h = np.zeros_like(mx, dtype=np.float32)
    mask = df > 0
    # Avoid division by zero; compute hue sectors
    r_is_max = (mx == r_f) & mask
    g_is_max = (mx == g_f) & mask
    b_is_max = (mx == b_f) & mask

    # Note: vectorized hue (0..6) then /6
    h[r_is_max] = ( (g_f[r_is_max] - b_f[r_is_max]) / df[r_is_max] ) % 6.0
    h[g_is_max] = ( (b_f[g_is_max] - r_f[g_is_max]) / df[g_is_max] ) + 2.0
    h[b_is_max] = ( (r_f[b_is_max] - g_f[b_is_max]) / df[b_is_max] ) + 4.0
    h *= (1.0 / 6.0)

    order_ring_colors_local = np.argsort(h, kind='mergesort')
    ring_color_ids = part_idx[order_ring_colors_local]
else:
    ring_color_ids = np.array([], dtype=np.uint32)

# The remaining colors go to non-ring pixels, ordered by luminance (dark -> bright).
# Integer luma approximation: Y ≈ 0.2126R + 0.7152G + 0.0722B, scaled by 10,000 for integers.
if N_non > 0:
    used_mask = np.zeros(N, dtype=bool)
    used_mask[ring_color_ids] = True
    non_color_ids = np.flatnonzero(~used_mask)

    # Luma as integer to avoid huge float arrays
    Yint = (2126 * R[non_color_ids].astype(np.uint32) +
            7152 * G[non_color_ids].astype(np.uint32) +
             722 * B[non_color_ids].astype(np.uint32))  # 0..2,550,000

    order_non_local = np.argsort(Yint, kind='mergesort')
    nonring_color_ids = non_color_ids[order_non_local]
else:
    nonring_color_ids = np.array([], dtype=np.uint32)

# --- 3) Assemble the final image buffer (bijective mapping) ---

img = np.empty((N_total, 3), dtype=np.uint8)

# Assign ring: pixels ordered by azimuth, colors ordered by hue
if N_ring > 0:
    Rc = (ring_color_ids >> 16) & 0xFF
    Gc = (ring_color_ids >>  8) & 0xFF
    Bc = (ring_color_ids      ) & 0xFF
    img[order_ring_pixels, 0] = Rc.astype(np.uint8)
    img[order_ring_pixels, 1] = Gc.astype(np.uint8)
    img[order_ring_pixels, 2] = Bc.astype(np.uint8)

# Assign non-ring: pixels ordered by intensity (dark->bright), colors ordered by luma
if N_non > 0:
    Rn = (nonring_color_ids >> 16) & 0xFF
    Gn = (nonring_color_ids >>  8) & 0xFF
    Bn = (nonring_color_ids      ) & 0xFF
    img[order_non_pixels, 0] = Rn.astype(np.uint8)
    img[order_non_pixels, 1] = Gn.astype(np.uint8)
    img[order_non_pixels, 2] = Bn.astype(np.uint8)

# Reshape and save
img = img.reshape((H, W, 3))
Image.fromarray(img, mode="RGB").save(OUTFILE, compress_level=0)

# --- 4) Verify exact usage of all 16,777,216 colors ---
packed = (img[...,0].astype(np.uint32) << 16) | (img[...,1].astype(np.uint32) << 8) | img[...,2].astype(np.uint32)
uniq = np.unique(packed)
assert uniq.size == (1 << 24), f"Color count mismatch: {uniq.size}"
print("Success!")
print(f"Ring pixels: {N_ring:,} | Non-ring pixels: {N_non:,}")
print(f"Saved: {OUTFILE}")

Why this maps to the real M87*

• The ring with a southern bright sector and deep central depression mirrors what EHT reported (ring diameter ≈ 42 μas, strong azimuthal asymmetry consistent with relativistic beaming). The 2023 PRIMO reconstruction tightened the ring width; we emulate a thin band. The 2024 multi-epoch analysis again found a southerly brightness distribution and a ring size consistent with 42 μas.