Miura 折叠¶

"""三浦折叠：一拉即折拢整张薄片。

一张折成相同平行四边形（边 a、b，角 γ）的薄片可刚性折叠——无任何面片弯曲或拉伸——仅
有单一自由度，即折叠角 θ。采用 Schenk 与 Guest 的单元胞（PNAS 2013）：
H = a sin θ sin γ，S = b cos θ tan γ / √(1 + cos²θ tan²γ)、
L = a √(1 − sin²θ sin²γ)，V = b / √(1 + cos²θ tan²γ)，顶点 (i, j) 位于
(i S, j L + (i mod 2) V, (j mod 2) H)：每条边保持其长度，每个面片保持平面（在构建时
检验）。宽度与长度同时收缩，泊松比为负 ν = −tan²γ cos²θ，这正是三浦公夫的折叠能以一拉
展开卫星太阳能板的原因（Space Flyer Unit，1995）。
"""

import numpy as np

import manimgx as m


def camera_axes(phi: float, theta: float, gamma: float) -> np.ndarray:
    """相机轨道角对应的坐标轴，以行向量形式：右、上、朝向观察者。"""

    def spin(a: float) -> np.ndarray:  # 绕 z 轴旋转
        return np.array(
            [[np.cos(a), -np.sin(a), 0.0], [np.sin(a), np.cos(a), 0.0], [0.0, 0.0, 1.0]]
        )

    c, s = np.cos(phi), np.sin(phi)
    tilt = np.array(
        [[1.0, 0.0, 0.0], [0.0, c, s], [0.0, -s, c]]
    )  # 绕 x 轴旋转 −φ
    return spin(gamma) @ tilt @ spin(-theta - np.pi / 2)


A, B, GAMMA = 1.0, 1.0, np.radians(60)  # 面片边长与夹角
COLS, ROWS = 14, 12  # 横向面片数（b 边）与纵向面片数（a 边）
SCALE = 0.44  # 薄片单位 → 场景单位
PAPER, INK = "#d3c8ae", "#2f4b7c"  # 薄片颜色，折叠后交替列的墨色
MOUNTAIN, VALLEY = "#d62828", "#1d5fd8"


def cell(theta: float) -> tuple[float, float, float, float]:
    """折叠角 θ 下 Schenk 与 Guest 的单元胞：高度 H、锯齿步距 S、行间距 L、偏移 V。"""
    s, c = np.sin(theta), np.cos(theta)
    k = np.sqrt(1 + (c * np.tan(GAMMA)) ** 2)
    return (
        A * s * np.sin(GAMMA),
        B * c * np.tan(GAMMA) / k,
        A * np.sqrt(1 - (s * np.sin(GAMMA)) ** 2),
        B / k,
    )


def sheet(theta: float) -> np.ndarray:
    """折叠后薄片的顶点 (COLS + 1, ROWS + 1, 3)，以薄片单位计。"""
    h, s, pitch, v = cell(theta)
    i, j = np.meshgrid(np.arange(COLS + 1), np.arange(ROWS + 1), indexing="ij")
    return np.stack([i * s, j * pitch + (i % 2) * v, (j % 2) * h], -1).astype(float)


def rigidity(theta: float) -> float:
    """折叠角 θ 下，任一边长或任一面片平面度的最大误差。"""
    p = sheet(theta)
    across = np.linalg.norm(p[1:] - p[:-1], axis=-1) - B
    along = np.linalg.norm(p[:, 1:] - p[:, :-1], axis=-1) - A
    o, x, y, far = p[:-1, :-1], p[1:, :-1], p[:-1, 1:], p[1:, 1:]
    normal = np.cross(x - o, y - o)
    normal /= np.linalg.norm(normal, axis=-1, keepdims=True)
    warp = np.einsum("ijk,ijk->ij", far - o, normal)
    return float(max(np.abs(across).max(), np.abs(along).max(), np.abs(warp).max()))


def facets(p: np.ndarray) -> np.ndarray:
    """每个面片各自的四个角点 (COLS · ROWS · 4, 3)：不共享，使每个面片平面着色。"""
    quads = np.stack([p[:-1, :-1], p[1:, :-1], p[1:, 1:], p[:-1, 1:]], 2)
    return quads.reshape(-1, 3)


def creases() -> tuple[np.ndarray, np.ndarray]:
    """每条内部折痕以 (顶点 a, 顶点 b, 面片 1, 面片 2) 索引行指向顶点网格和面片列表，
    并标记是否为山形折痕（从上方看）。"""
    vid = np.arange((COLS + 1) * (ROWS + 1)).reshape(COLS + 1, ROWS + 1)
    fid = np.arange(COLS * ROWS).reshape(COLS, ROWS)
    rows = []
    for i in range(COLS):  # 面片行 j − 1 与 j 之间的锯齿折痕
        for j in range(1, ROWS):
            rows.append([vid[i, j], vid[i + 1, j], fid[i, j - 1], fid[i, j]])
    for i in range(1, COLS):  # 面片列 i − 1 与 i 之间的折痕
        for j in range(ROWS):
            rows.append([vid[i, j], vid[i, j + 1], fid[i - 1, j], fid[i, j]])
    index = np.array(rows)
    # 折痕为山形当其高于两侧面片中心连线的中点
    p = sheet(np.radians(10)).reshape(-1, 3)
    centers = facets(sheet(np.radians(10))).reshape(-1, 4, 3).mean(axis=1)
    middle = (p[index[:, 0]] + p[index[:, 1]]) / 2
    return index, middle[:, 2] > (centers[index[:, 2], 2] + centers[index[:, 3], 2]) / 2


def ribbons(p: np.ndarray, index: np.ndarray, width: np.ndarray) -> np.ndarray:
    """沿各折痕宽 `width` 的条带（每条折痕一条），在两侧面片上各占一半，
    略微抬离面片：(creases · 8, 3) 个顶点，每半 4 个。"""
    grid = p.reshape(-1, 3)
    corners = facets(p).reshape(-1, 4, 3)
    a, b = grid[index[:, 0]], grid[index[:, 1]]
    halves = []
    for side in (2, 3):
        center = corners[index[:, side]].mean(axis=1)
        normal = np.cross(
            corners[index[:, side], 1] - corners[index[:, side], 0],
            corners[index[:, side], 3] - corners[index[:, side], 0],
        )
        normal /= np.linalg.norm(normal, axis=1, keepdims=True)
        along = (b - a) / np.linalg.norm(b - a, axis=1, keepdims=True)
        inward = center - a - ((center - a) * along).sum(axis=1, keepdims=True) * along
        inward /= np.linalg.norm(inward, axis=1, keepdims=True)
        lift = 0.012 * normal
        side_edge = width[:, None] * inward
        halves.append(np.stack([a, b, b + side_edge, a + side_edge], 1) + lift[:, None])
    return np.concatenate(halves, 1).reshape(-1, 3)


def project(camera: m.Camera, point: np.ndarray) -> np.ndarray:
    """三维点在屏幕上的落点，以帧坐标表示（用于平视显示标签）。"""
    turned = camera_axes(camera.get_phi(), camera.get_theta(), camera.get_gamma())
    q = turned @ (point - camera.frame_center)
    depth = 1 - q[2] / camera.get_focal_distance()
    return np.array([*(camera.get_zoom() * q[:2] / depth), 0.0])


def rgba(color: str) -> np.ndarray:
    return np.array([*m.ManimColor(color).to_rgb(), 1.0])


class MiuraOri(m.ThreeDScene):
    def construct(self) -> None:
        self.camera.background_color = "#0b1020"
        self.camera.light_source.move_to([-6, -9, 7])
        error = max(rigidity(np.radians(d)) for d in np.linspace(0, 89, 90))
        fold = m.ValueTracker(0.0)  # 唯一自由度，以度计
        crease_width = m.ValueTracker(0.03)
        sweep = m.ValueTracker(0.0)  # 折痕沿薄片逐渐绘出

        def now(extra: float = 0.0) -> np.ndarray:
            """当前折叠角下（加 `extra` 度）的薄片，以场景单位计，居中。"""
            p = sheet(np.radians(fold.get_value() + extra))
            middle = (p.reshape(-1, 3).max(axis=0) + p.reshape(-1, 3).min(axis=0)) / 2
            return SCALE * (p - middle)

        # 薄片，每个面片平面着色；折叠时交替列渐变为墨色，使两种
        # 面片倾斜方向得以显现
        odd = np.repeat(np.arange(COLS * ROWS) // ROWS % 2 == 1, 4)[:, None]
        quad = np.array([[0, 1, 2], [0, 2, 3]])
        paper = m.MeshMobject(
            facets(now()),
            (np.arange(COLS * ROWS)[:, None, None] * 4 + quad).reshape(-1, 3),
            shade_in_3d=True,
        )

        def refold(mob: m.MeshMobject) -> None:
            mob.points = facets(now())
            ink = np.clip((fold.get_value() - 8) / 30, 0, 1)
            ink = ink * ink * (3 - 2 * ink)
            tint = rgba(PAPER) + ink * (rgba(INK) - rgba(PAPER))
            mob.paint = mob.paint.but(fill=np.where(odd, tint, rgba(PAPER)))

        refold(paper)
        paper.add_updater(refold)

        # 山形与谷形折痕，绘为随折叠加深而变细的条带
        index, mountain = creases()
        colors = np.where(np.repeat(mountain, 2)[:, None], rgba(MOUNTAIN), rgba(VALLEY))
        flat = now().reshape(-1, 3)
        across = (flat[index[:, 0], :2] + flat[index[:, 1], :2]).sum(axis=1)
        across = (across - across.min()) / (across.max() - across.min())

        def widths() -> np.ndarray:
            drawn = np.clip((1.3 * sweep.get_value() - across) / 0.3, 0, 1)
            return crease_width.get_value() * drawn

        lines = m.MeshMobject(
            ribbons(now(), index, widths()),
            (np.arange(len(index) * 2)[:, None, None] * 4 + quad).reshape(-1, 3),
            vertex_colors=np.repeat(colors, 4, axis=0),
            shade_in_3d=True,
        )
        lines.add_updater(
            lambda mob: setattr(mob, "points", ribbons(now(), index, widths()))
        )

        # 相机将薄片中心保持在屏幕上某固定点，并随薄片收缩
        # 而放大（按其测得的对角线）
        spot = np.array([2.5, -0.3])

        def span() -> float:
            p = now()
            return float(np.hypot(p[-1, 0, 0] - p[0, 0, 0], p[0, -1, 1] - p[0, 0, 1]))

        flat_span = span()

        def follow(mob: m.Mobject) -> None:
            cam = self.camera
            cam.zoom_tracker.set_value(0.95 * (flat_span / span()) ** 0.85)
            turned = camera_axes(cam.get_phi(), cam.get_theta(), cam.get_gamma())
            offset = (spot[0] * turned[0] + spot[1] * turned[1]) / cam.get_zoom()
            cam.frame.move_to(-offset)

        # 从所绘薄片测量：沿列方向的宽度、沿行方向的长度，
        # 泊松比由再多折一点的两状态求得
        def measure(extra: float = 0.0) -> tuple[float, float]:
            p = now(extra)
            return float(p[-1, 0, 0] - p[0, 0, 0]), float(p[0, -1, 1] - p[0, 0, 1])

        flat_w, flat_l = measure()  # 薄片初始为平面

        def poisson() -> float:
            (w0, l0), (w1, l1) = measure(), measure(0.01)
            return -np.log(l1 / l0) / np.log(w1 / w0)

        # 拉力：沿宽度向外（蓝色）；薄片沿长度方向响应（金色）
        pulling = m.ValueTracker(0.0)

        def arrows() -> m.VGroup:
            p = now().reshape(-1, 3)
            low, high = p.min(axis=0), p.max(axis=0)
            mid = (low + high) / 2
            out = m.VGroup()
            for axis, color in ((0, "#8ecae6"), (1, "#ffb703")):
                for side, edge in ((-1, low), (1, high)):
                    start = mid.copy()
                    start[axis] = edge[axis] + side * 0.15
                    stop = start.copy()
                    stop[axis] += side * 0.9
                    a, b = project(self.camera, start), project(self.camera, stop)
                    out.add(m.Arrow(a, b, buff=0, stroke_width=6, color=color))
            return out.set_opacity(pulling.get_value())

        pull = m.always_redraw(arrows)

        # 平视显示
        title = m.Text("Miura-ori", font_size=40).to_corner(m.UL)
        notes = [
            "a sheet creased into identical parallelograms",
            "one angle folds every crease at once, rigidly",
            "pull it apart one way: it opens the other way too",
        ]
        subtitles = [m.Text(t, font_size=24).set_color(m.GREY_B) for t in notes]
        for sub in subtitles:
            sub.next_to(title, m.DOWN, aligned_edge=m.LEFT, buff=0.15)
        legend = m.VGroup()
        for name, color in (("mountain fold", MOUNTAIN), ("valley fold", VALLEY)):
            swatch = m.Line(m.ORIGIN, 0.5 * m.RIGHT, stroke_width=6, color=color)
            legend.add(
                m.VGroup(swatch, m.Text(name, font_size=22)).arrange(m.RIGHT, buff=0.2)
            )
        legend.arrange(m.DOWN, aligned_edge=m.LEFT, buff=0.2)
        legend.move_to([-6.7, 2.15, 0], aligned_edge=m.LEFT)
        shown = m.ValueTracker(0.0)

        def dial() -> m.VGroup:
            """折叠角以滑块表示，宽度与长度以条形显示，均为实测值。"""
            fade = shown.get_value()
            track = m.Line(
                [-6.7, 1.2, 0], [-3.7, 1.2, 0], stroke_width=3, color=m.GREY_C
            )
            at = -6.7 + 3.0 * fold.get_value() / 90
            knob = m.Dot([at, 1.2, 0], 0.09, color=m.WHITE)
            angle = m.MathTex(rf"\theta = {fold.get_value():.0f}^\circ", font_size=30)
            angle.next_to(track, m.RIGHT, 0.3)
            width, length = measure()
            bars = m.VGroup()
            for k, (name, ratio, color) in enumerate(
                (
                    ("width", width / flat_w, "#8ecae6"),
                    ("length", length / flat_l, "#ffb703"),
                )
            ):
                y = 0.45 - 0.55 * k
                word = m.Text(name, font_size=22).move_to(
                    [-6.7, y, 0], aligned_edge=m.LEFT
                )
                box = m.Rectangle(width=2.2 * ratio, height=0.22).set_fill(color, 0.9)
                box.set_stroke(width=0).move_to([-5.45, y, 0], aligned_edge=m.LEFT)
                value = m.Text(f"{100 * ratio:.0f}%", font_size=22).next_to(
                    box, m.RIGHT, 0.15
                )
                bars.add(word, box, value)
            nu = m.MathTex(
                rf"\nu = -\frac{{d L / L}}{{d W / W}} = {poisson():.2f}", font_size=32
            ).move_to([-6.7, -0.95, 0], aligned_edge=m.LEFT)
            group = m.VGroup(track, knob, angle, bars, nu)
            return group.set_opacity(fade)

        panel = m.always_redraw(dial)
        law = m.MathTex(r"\nu = -\tan^2\gamma\,\cos^2\theta", font_size=32)
        law.move_to([-6.7, -1.8, 0], aligned_edge=m.LEFT)
        mantissa, power = f"{error:.0e}".split("e")
        note = m.MathTex(
            rf"\text{{rigid: edges and flatness exact to }} {mantissa} \times"
            rf" 10^{{{int(power)}}}",
            font_size=24,
        ).set_color(m.GREY_B)
        note.move_to([-6.7, -2.45, 0], aligned_edge=m.LEFT)
        closing = m.VGroup(
            m.Text("Koryo Miura's fold opens", font_size=28),
            m.Text("a solar array with one pull", font_size=28),
        ).arrange(m.DOWN, aligned_edge=m.LEFT, buff=0.12)
        closing.move_to([-6.7, -3.3, 0], aligned_edge=m.LEFT)
        self.add_fixed_in_frame_mobjects(
            title, *subtitles, legend, panel, law, note, closing, pull
        )
        self.remove(*subtitles, legend, law, note, closing)

        driver = m.Mobject().add_updater(follow)
        self.set_camera_orientation(phi=54 * m.DEGREES, theta=-104 * m.DEGREES)
        self.begin_ambient_camera_rotation(rate=0.03)
        follow(driver)
        self.add(paper, lines, driver)
        # 0–4.5 秒：平面薄片及其折痕
        self.play(m.FadeIn(subtitles[0]), run_time=1.0)
        self.play(sweep.animate.set_value(1.0), m.FadeIn(legend), run_time=1.8)
        self.play(shown.animate.set_value(1.0), run_time=1.0)
        self.wait(0.8)
        # 4.5–14 秒：折叠，所有折痕同时动作
        self.play(m.FadeOut(subtitles[0]), run_time=0.3)
        self.play(m.FadeIn(subtitles[1]), run_time=0.4)
        self.play(
            fold.animate.set_value(82),
            crease_width.animate(rate_func=lambda t: min(1.0, 3 * t)).set_value(0.0),
            m.FadeOut(legend),
            run_time=8.5,
        )
        self.wait(1.0)
        # 15–23 秒：沿宽度拉开，长度方向也随之增长
        self.play(m.FadeOut(subtitles[1]), run_time=0.3)
        self.play(m.FadeIn(subtitles[2]), pulling.animate.set_value(1.0), run_time=0.6)
        self.play(fold.animate.set_value(28), m.FadeIn(law), run_time=6.5)
        self.play(pulling.animate.set_value(0.0), run_time=0.6)
        # 23–30.5 秒：定版画面
        self.play(
            m.FadeIn(note), m.FadeIn(closing), fold.animate.set_value(45), run_time=3
        )
        self.wait(30.5 - self.time)


if __name__ == "__main__":
    MiuraOri().render("miura_ori.mp4")