sciencehub-team/black-holePublic

# /// script
# requires-python = ">=3.12"
# dependencies = ["marimo>=0.25.1", "numpy", "matplotlib"]
# ///

import marimo

__generated_with = "0.25.1"
app = marimo.App(width="medium")


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        # What a black hole looks like

        Light passing a black hole follows the geometry of spacetime around it.
        For a non-rotating black hole of mass $M$ that geometry is the
        Schwarzschild solution, and a light ray in its equatorial plane obeys one
        equation (units with $G = c = 1$, $u = M/r$, angle $\varphi$):

        $$\frac{d^2u}{d\varphi^2} + u = 3u^2$$

        Without the $3u^2$ term this is a straight line. With it, rays bend, some
        circle the hole before escaping, and every ray aimed closer than a
        critical distance falls in. This notebook solves that equation and
        uses it twice: to draw what a camera near the hole would see, and to
        follow single rays.
        """
    )
    return


@app.cell
def _():
    import marimo as mo
    import numpy as np
    import matplotlib.pyplot as plt
    return mo, np, plt


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        ## The view from nearby

        Surround the hole with a thin, glowing disc of gas, the way matter
        spirals into real black holes, and trace one ray backwards from each
        pixel of a camera. Every ray lies in a plane through the hole, so the
        same equation decides where it meets the disc. Rays that miss the disc
        and escape show the sky behind, displaced by the same bending.

        The left side is brighter because the gas there moves towards the
        camera (Doppler beaming); the arch over the top and the ring
        underneath are the far side of the disc, seen around the hole. This is
        the calculation J.-P. Luminet made by hand in 1979.
        """
    )
    return


@app.cell
def _(np):
    def ray_table(b_max=30.0, nb=1500, step=0.002, max_angle=3.2 * np.pi):
        """r(phi) for rays arriving from infinity with impact parameter b (units of M).
        Returns b grid, phi grid, and u = M/r (nan once the ray is inside the horizon)."""
        b = np.linspace(0.2, b_max, nb)
        n = int(max_angle / step) + 1
        u = np.zeros(nb); du = 1 / b
        out = np.full((n, nb), np.nan, dtype=np.float32)
        out[0] = 0
        alive = np.ones(nb, bool)
        escape = np.full(nb, np.nan)
        prev = u.copy()
        for i in range(1, n):
            f = lambda x: 3 * x * x - x
            k1u, k1v = du, f(u)
            k2u, k2v = du + step/2*k1v, f(u + step/2*k1u)
            k3u, k3v = du + step/2*k2v, f(u + step/2*k2u)
            k4u, k4v = du + step*k3v, f(u + step*k3u)
            u = np.where(alive, u + step/6*(k1u+2*k2u+2*k3u+k4u), 0.5)
            du = np.where(alive, du + step/6*(k1v+2*k2v+2*k3v+k4v), 0.0)
            crossed = alive & np.isnan(escape) & (u < 0) & (i > 1)
            escape[crossed] = (i - 1 + prev[crossed] / (prev[crossed] - u[crossed])) * step
            alive &= (u < 0.5) & (u >= 0)
            prev = u.copy()
            out[i] = np.where(alive, u, np.nan)
        return b, step, out, escape

    def render(width, height, incl_deg=80, fov=34.0, r_in=6.0, r_out=24.0, table=None, ss=1, stars=None):
        b_grid, step, U, escape = table
        W, H = width*ss, height*ss
        x = (np.arange(W) - W/2 + 0.5) / W * fov
        y = (np.arange(H) - H/2 + 0.5) / W * fov
        X, Y = np.meshgrid(x, -y)
        b = np.hypot(X, Y); th = np.arctan2(Y, X)
        i = np.radians(incl_deg)
        phi0 = np.mod(np.arctan2(-np.cos(i), np.sin(i)*np.sin(th)), np.pi)
        img = np.zeros_like(b)
        bi = (b - b_grid[0]) / (b_grid[1]-b_grid[0])
        valid_b = (bi >= 0) & (bi < len(b_grid)-1)
        bi = np.clip(bi, 0, len(b_grid)-1.001); b0 = bi.astype(int); fb = bi - b0
        for k in (2, 1, 0):  # far images first, nearest image drawn last
            phi = phi0 + k*np.pi
            pi_ = phi / step
            ok = valid_b & (pi_ < U.shape[0]-1)
            pi_ = np.clip(pi_, 0, U.shape[0]-1.001); p0 = pi_.astype(int); fp = pi_ - p0
            u = (U[p0, b0]*(1-fp)*(1-fb) + U[p0+1, b0]*fp*(1-fb) + U[p0, b0+1]*(1-fp)*fb + U[p0+1, b0+1]*fp*fb)
            with np.errstate(invalid='ignore', divide='ignore'):
                r = 1/u
                hit = ok & np.isfinite(r) & (r >= r_in) & (r <= r_out)
                # Keplerian emitter; photon L_z = -alpha sin i, alpha = b cos(theta)
                omega = r**-1.5
                g = np.sqrt(1 - 3/r) / (1 + omega * X * np.sin(i))
                flux = (1 - np.sqrt(r_in/r)) / r**3   # Page-Thorne-like radial profile
                I = flux * g**4
            img = np.where(hit, I, img)
        # Sky: rays that miss the disc and escape leave along the swept angle phi_f.
        sky = np.zeros_like(img)
        if stars is not None:
            phi_f = np.interp(b, b_grid, np.nan_to_num(escape, nan=-1.0))
            free = (img == 0) & (phi_f > 0) & (b < b_grid[-1])
            n_obs = np.array([0.0, -np.sin(i), np.cos(i)])  # observer direction, disc normal = z
            e_x = np.array([1.0, 0.0, 0.0])
            e_y = np.array([0.0, np.cos(i), np.sin(i)])
            ex = np.cos(th)[..., None] * e_x + np.sin(th)[..., None] * e_y
            d = -(np.cos(phi_f)[..., None] * n_obs + np.sin(phi_f)[..., None] * ex)
            lon = np.arctan2(d[..., 1], d[..., 0]); lat = np.arcsin(np.clip(d[..., 2], -1, 1))
            h, w = stars.shape
            px = ((lon + np.pi) / (2 * np.pi) * w).astype(int) % w
            py = np.clip(((np.pi / 2 - lat) / np.pi * h).astype(int), 0, h - 1)
            sky = np.where(free, stars[py, px], 0.0)
        if ss > 1:
            img = img.reshape(height, ss, width, ss).mean(axis=(1, 3))
            sky = sky.reshape(height, ss, width, ss).mean(axis=(1, 3))
        return img, sky

    def star_texture(w=4096, h=2048, n=9000, seed=7):
        """Equirectangular sky with stars spread uniformly over the sphere."""
        rng = np.random.default_rng(seed)
        z = rng.uniform(-1, 1, n); lon = rng.uniform(-np.pi, np.pi, n)
        lat = np.arcsin(z)
        x = ((lon + np.pi) / (2 * np.pi) * w).astype(int) % w
        y = np.clip(((np.pi / 2 - lat) / np.pi * h).astype(int), 0, h - 1)
        tex = np.zeros((h, w), np.float32)
        np.add.at(tex, (y, x), rng.pareto(1.6, n).astype(np.float32) * 0.15 + 0.05)
        for dy, dx in ((0, 1), (1, 0)):
            tex += 0.35 * np.roll(np.roll(tex, dy, 0), dx, 1)
        return np.clip(tex, 0, 1)

    RAY_TABLE = ray_table(b_max=60, nb=1200, step=0.003, max_angle=3.3 * np.pi)
    STARS = star_texture()
    return RAY_TABLE, STARS, ray_table, render, star_texture


@app.cell
def _(mo):
    inclination = mo.ui.slider(0, 89, step=1, value=84, label="Camera angle above the disc axis (degrees)")
    inclination
    return (inclination,)


@app.cell
def _(RAY_TABLE, STARS, inclination, mo, np, plt, render):
    import io
    from matplotlib import cm

    disc, sky = render(640, 360, incl_deg=inclination.value, fov=44, r_out=18, table=RAY_TABLE, ss=1, stars=STARS)
    level = np.clip(disc / np.percentile(disc[disc > 0], 99.7), 0, 1) ** 0.55
    picture = np.clip(cm.afmhot(level / 1.1)[..., :3] + np.clip(sky, 0, 1)[..., None] * np.array([0.85, 0.9, 1.0]), 0, 1)
    buffer = io.BytesIO()
    plt.imsave(buffer, picture, format="png")
    mo.image(buffer.getvalue(), width="100%")
    return


@app.cell
def _(np):
    B_CRITICAL = 3 * np.sqrt(3)  # impact parameter of the photon sphere, in units of M
    R_START = 60.0  # rays start this far from the hole, in units of M

    def trace(b, step=0.004, max_angle=6 * np.pi):
        """Trace rays with impact parameters `b` and return their paths.

        Integrates u'' = 3u^2 - u with fourth-order Runge-Kutta in the angle
        swept since the start, all rays at once. A ray ends when it crosses the
        horizon (r = 2M) or leaves the far boundary.
        """
        b = np.atleast_1d(np.asarray(b, dtype=float))
        u = np.full_like(b, 1 / R_START)
        du = np.sqrt(np.maximum(1 / b**2 - u**2 * (1 - 2 * u), 0))
        psi = np.arcsin(np.minimum(b / R_START, 1.0))
        alive = np.ones_like(b, dtype=bool)
        captured = np.zeros_like(b, dtype=bool)
        history_u, history_psi, history_alive = [u], [psi], [alive]

        def accel(u):
            return 3 * u**2 - u

        while alive.any() and (psi[alive] < max_angle).any():
            k1u, k1v = du, accel(u)
            k2u, k2v = du + step / 2 * k1v, accel(u + step / 2 * k1u)
            k3u, k3v = du + step / 2 * k2v, accel(u + step / 2 * k2u)
            k4u, k4v = du + step * k3v, accel(u + step * k3u)
            u = np.where(alive, u + step / 6 * (k1u + 2 * k2u + 2 * k3u + k4u), u)
            du = np.where(alive, du + step / 6 * (k1v + 2 * k2v + 2 * k3v + k4v), du)
            psi = np.where(alive, psi + step, psi)
            fell = alive & (u >= 0.5)
            left = alive & (u <= 1 / R_START) & (du < 0)
            captured |= fell
            history_u.append(u)
            history_psi.append(psi)
            history_alive.append(alive)
            alive = alive & ~(fell | left) & (psi < max_angle)

        r = 1 / np.maximum(np.array(history_u), 1e-9)
        angle = np.array(history_psi)
        drawn = np.array(history_alive)
        paths = [
            np.column_stack([-r[drawn[:, i], i] * np.cos(angle[drawn[:, i], i]),
                             r[drawn[:, i], i] * np.sin(angle[drawn[:, i], i])])
            for i in range(len(b))
        ]
        return paths, captured, psi

    return B_CRITICAL, trace


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        """
        ## Single rays

        Move the slider to aim a ray past the hole.
        """
    )
    return


@app.cell
def _(B_CRITICAL, mo):
    aim = mo.ui.slider(
        start=2.0,
        stop=12.0,
        step=0.01,
        value=round(B_CRITICAL + 0.02, 2),
        label="Impact parameter b (in units of M)",
        full_width=True,
    )
    aim
    return (aim,)


@app.cell
def _(np, trace):
    fan_paths, fan_captured, _ = trace(np.linspace(1.5, 11.5, 24))
    return fan_captured, fan_paths


@app.cell
def _(aim, fan_captured, fan_paths, np, plt, trace):
    (ray,), (ray_captured,), (ray_angle,) = [x[:1] for x in trace([aim.value])]

    fig, ax = plt.subplots(figsize=(7, 5.2))
    for path, fell in zip(fan_paths, fan_captured):
        ax.plot(*path.T, lw=0.7, color="#9a6700" if fell else "#8c959f", alpha=0.55)
    ax.plot(*ray.T, lw=2.2, color="#cf222e" if ray_captured else "#0969da")
    ax.add_patch(plt.Circle((0, 0), 2, color="black", zorder=3))
    ax.add_patch(plt.Circle((0, 0), 3, fill=False, ls="--", color="#6e7781", lw=0.8))
    ax.set(xlim=(-22, 22), ylim=(-14, 16), aspect="equal", xlabel="x / M", ylabel="y / M")
    ax.set_title(
        f"b = {aim.value:.2f} M: "
        + (
            "falls in"
            if ray_captured
            else f"escapes after turning {np.degrees(ray_angle - np.pi):.0f}° away from straight"
        )
    )
    fig.tight_layout()
    fig
    return


@app.cell(hide_code=True)
def _(B_CRITICAL, mo):
    mo.md(
        f"""
        The black disc is the event horizon at $r = 2M$. The dashed circle is the
        photon sphere at $r = 3M$, where light can orbit the hole. Rays aimed with
        $b < 3\\sqrt{{3}}\\,M \\approx {B_CRITICAL:.3f}\\,M$ (amber) fall in; rays
        aimed just above it wind around the photon sphere before escaping. Set
        the slider to 5.20 and then to 5.19 to see the boundary.
        """
    )
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        ## Deflection angle against Einstein's prediction

        Far from the hole the deflection is small, and Einstein's weak-field
        result $\delta = 4M/b$ holds. That is the figure Eddington's 1919
        eclipse expedition tested for light grazing the Sun. Close to the hole
        the full equation departs from it and the deflection grows without
        bound as $b$ approaches the critical value.
        """
    )
    return


@app.cell
def _(B_CRITICAL, np, plt):
    def deflection(b, step=0.002):
        """Total bending of rays that come from and return to infinity (u = 0)."""
        u = np.zeros_like(b)
        du = 1 / b
        psi = np.zeros_like(b)
        out = np.full_like(b, np.nan)
        while np.isnan(out).any() and psi[0] < 8 * np.pi:
            k1u, k1v = du, 3 * u**2 - u
            k2u, k2v = du + step / 2 * k1v, 3 * (u + step / 2 * k1u) ** 2 - (u + step / 2 * k1u)
            k3u, k3v = du + step / 2 * k2v, 3 * (u + step / 2 * k2u) ** 2 - (u + step / 2 * k2u)
            k4u, k4v = du + step * k3v, 3 * (u + step * k3u) ** 2 - (u + step * k3u)
            new_u = u + step / 6 * (k1u + 2 * k2u + 2 * k3u + k4u)
            du = du + step / 6 * (k1v + 2 * k2v + 2 * k3v + k4v)
            crossed = np.isnan(out) & (new_u < 0) & (psi > 0)
            out[crossed] = psi[crossed] + step * u[crossed] / (u[crossed] - new_u[crossed])
            u, psi = new_u, psi + step
        return out - np.pi

    impact = B_CRITICAL + np.geomspace(0.005, 95, 50)
    bent = np.degrees(deflection(impact))

    fig2, ax2 = plt.subplots(figsize=(7, 3.6))
    ax2.loglog(impact, bent, color="#0969da", lw=2, label="Schwarzschild, integrated")
    ax2.loglog(impact, np.degrees(4 / impact), color="#8c959f", ls="--", label="4M/b, weak field")
    ax2.axvline(B_CRITICAL, color="#9a6700", lw=0.8)
    ax2.set(xlabel="impact parameter b / M", ylabel="deflection (degrees)")
    ax2.legend(frameon=False)
    fig2.tight_layout()
    fig2
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        ## The shadow of M87*

        A distant observer sees a dark disc where rays with $b < 3\sqrt{3}\,M$
        would have come from. Its angular diameter is
        $2 \cdot 3\sqrt{3}\, GM / (c^2 D)$. For the black hole in the galaxy
        M87, with a mass of about 6.5 billion Suns at 16.8 megaparsecs (the
        values the Event Horizon Telescope used in 2019), that gives:
        """
    )
    return


@app.cell
def _(B_CRITICAL, mo, np):
    G, c = 6.674e-11, 2.998e8
    M_SUN, MPC = 1.989e30, 3.086e22
    MICROARCSEC = np.radians(1 / 3.6e9)

    mass = mo.ui.number(start=0.1, stop=20, step=0.1, value=6.5, label="Mass (billion Suns)")
    distance = mo.ui.number(start=1, stop=100, step=0.1, value=16.8, label="Distance (Mpc)")
    mo.hstack([mass, distance], justify="start")
    return G, MICROARCSEC, MPC, M_SUN, c, distance, mass


@app.cell
def _(B_CRITICAL, G, MICROARCSEC, MPC, M_SUN, c, distance, mass, mo):
    theta = 2 * B_CRITICAL * G * mass.value * 1e9 * M_SUN / (c**2 * distance.value * MPC)
    mo.md(
        f"""
        **Shadow diameter: {theta / MICROARCSEC:.1f} microarcseconds.** The Event
        Horizon Telescope measured a ring of 42 ± 3 microarcseconds
        (EHT Collaboration 2019, ApJL 875, L1). The ring sits slightly outside
        the shadow because the light we see is emitted near it.
        """
    )
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        """
        ---
        **Method.** Fourth-order Runge-Kutta in the swept angle, step 0.004 rad
        (0.002 for the deflection curve). The pictured rays start at $r = 60M$;
        the deflection curve integrates from and back to infinity ($u = 0$). **References.** K. Schwarzschild (1916); S. Chandrasekhar,
        *The Mathematical Theory of Black Holes* (1983), ch. 3; Event Horizon
        Telescope Collaboration, ApJL 875 L1 (2019); J.-P. Luminet, A&A 75 (1979) 228.
        The disc is thin and Keplerian, emits from $r = 6M$ to $18M$ with a
        Page-Thorne-like profile, and the star field is synthetic.
        Fork it to trace a rotating (Kerr) black hole or a thick, turbulent disc.
        """
    )
    return


if __name__ == "__main__":
    app.run()