sciencehub-team/black-holePublic

import marimo

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


@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, render


@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 _(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()