シミュレーションと３次元コンピュータグラフィックス

【サイト内の関連ページ】

• Panda3D 3次元ゲームエンジン入門: 別ページ »で説明
• Panda3Dゲームエンジン：システム構築の基礎 別ページ »で説明
• WindowsでのPanda3Dのインストール: 別ページ »で説明
• UbuntuでのPanda3Dのインストール: 別ページ »で説明
• Panda3Dの機能概要(説明資料)[PDF]、[パワーポイント]

共通の準備（Panda3Dのインストールとプログラムの実行）

Panda3Dのインストール

3つの演習はすべてPanda3Dを使う。Panda3Dは3次元コンピュータグラフィックスのためのライブラリである。インストールは1回だけ行えばよい。演習３はnumbaを使う。

手順：コマンドプロンプトを管理者権限で起動する（Windowsキー → 「cmd」と入力 → 右クリック → 「管理者として実行」）。起動したコマンドプロンプトに次を入力して実行する。

pip install --no-user -U panda3d numba

プログラムの実行

各演習のコードを実行する方法は2通りある。どちらでも動作する。

方法A：Visual Studio Codeで実行する。Visual Studio Codeがインストール・設定済みであれば、本ページのコードを編集画面にコピー＆ペーストし、そのまま実行する。

方法B：ファイルに保存して実行する。メモ帳などのテキストエディタを開き、本ページのコードを貼り付ける。a.pyのようなファイル名で保存する。コマンドプロンプトでそのファイルがあるフォルダに移動し、次を入力して実行する。

python a.py

共通のマウス操作

3つの演習はすべて同じマウス操作で視点を動かせる。

• 左ドラッグ：視点を回転
• 右ドラッグ：ズームイン・ズームアウト
• 中ドラッグ：視点を平行移動

3次元コンピュータグラフィックスの利点は、固定された1枚絵ではなく、見る側が自由に視点を変えて対象を観察できる点にある。各演習ではマウス操作で視点を変えながら観察することを推奨する。

重力N体シミュレーション（銀河衝突）

このプログラムが描くもの

2つの銀河（青25個・赤25個の天体）を左右から接近させて衝突させる。各天体は他のすべての天体から重力を受け、その合力に従って運動する。画面には次の3要素が描かれる。

• 点：各天体の現在位置（青または赤、6ピクセル）
• 線：過去500ステップ分の軌跡。古い部分ほど薄く（フェードアウト）
• 文字：画面左上に物理量（エネルギー、運動量）を表示

複雑な形状のモデルを用意しなくても、点と線だけで天文現象のような壮大な対象を可視化できる。本プログラムはその一例である。

背景となる考え方

ニュートンの万有引力の法則により、2つの物体間の引力は両者の質量の積に比例し、距離の2乗に反比例する。N個の天体があれば、すべてのペアについてこの計算を行い、各天体に働く力の合計から加速度を求め、速度と位置を時間刻みだけ進める。これを毎フレーム繰り返すことで、運動が時間発展する。

物理的には、外力が働かない孤立系では全エネルギーと全運動量が保存される。画面左上の表示でこれを確認できれば、シミュレーションが妥当に動いている証拠となる。

演習

《手順》

1. 共通の準備の手順でコードを実行する。2つの銀河が左右から接近する
2. 衝突前、衝突中、衝突後の3段階で、それぞれ30秒程度ずつ画面を観察する
3. マウス操作で視点を動かし、真上からと真横からの2方向でも観察する
4. 画面左上の「Total E」（全エネルギー）と「Momentum」（運動量の大きさ）を、実行開始直後・衝突中・衝突後の3時点で確認する

《ヒント》

• 衝突は実行開始からおよそ15秒～30秒後に起こる
• 軌跡の色（青／赤）で元の所属銀河が分かる。融合後にどちらの天体がどこへ行ったかを追える
• 真上から見ると円盤状の回転構造が分かりやすい。真横から見ると銀河の薄さ（z方向の広がりの小ささ）が分かる
• 「Total E」は負の値になる（位置エネルギーが負のため）。重要なのは絶対値ではなく、時間が経っても大きく変わらないこと

《考察ポイント》

• 軌跡フェードの効果：軌跡が古いものほど薄くなる表現は、運動の方向と速さを直感的に伝える。もし軌跡が一定の濃さだったら、画面はどう見えるか想像する
• 3次元で見ることの意味：2次元の図（教科書の挿絵）と、自分でマウス操作して視点を変えながら見る3次元表示では、衝突の理解にどんな違いがあるか
• 物理量の保存：Total EとMomentumがほぼ一定に保たれていたか。大きく変動した場合、数値計算の誤差が原因と考えられる
• ルールから形が生まれる：プログラムには「衝突後の銀河の形」を直接指定するコードは存在しない。それでも引力の法則だけで複雑な形状が生まれる。これは3次元コンピュータグラフィックスにおける「ルールが形を作る」設計思想の典型例である

# コード1：重力N体シミュレーション（銀河衝突）
from direct.showbase.ShowBase import ShowBase
from direct.task import Task
from panda3d.core import GeomVertexFormat, GeomVertexData, GeomVertexWriter
from panda3d.core import Mat4
from panda3d.core import Geom, GeomNode, GeomPoints, GeomLines
from panda3d.core import RenderModeAttrib, TextNode
from direct.gui.OnscreenText import OnscreenText
import numpy as np
import random
import math


class Body:
    """天体を表すクラス"""
    def __init__(self, x, y, z, vx, vy, vz, mass, color):
        self.pos = np.array([x, y, z], dtype=np.float64)
        self.vel = np.array([vx, vy, vz], dtype=np.float64)
        self.mass = mass
        self.color = color
        self.trail = []

    def add_trail_point(self):
        self.trail.append(self.pos.copy())
        if len(self.trail) > 500:
            self.trail.pop(0)


class GalaxyCollision(ShowBase):
    def __init__(self):
        ShowBase.__init__(self)
        self.setBackgroundColor(0, 0, 0)

        # 物理定数
        self.G = 1.0
        self.softening = 0.5
        self.dt = 0.01

        # 天体リスト
        self.bodies = []
        self.create_galaxies()

        # 天体メッシュ（点）
        self.particle_mesh = self.create_particle_mesh()
        node = self.render.attachNewNode(self.particle_mesh)
        node.setRenderMode(RenderModeAttrib.MPoint, 6.0)

        # 軌跡メッシュ（線）
        self.trail_mesh = self.create_trail_mesh()
        self.render.attachNewNode(self.trail_mesh)

        # カメラ設定
        self.disableMouse()
        self.camera.setPos(0, -100, 40)
        self.camera.lookAt(0, 0, 0)
        mat = Mat4(self.camera.getMat())
        mat.invertInPlace()
        self.mouseInterfaceNode.setMat(mat)
        self.enableMouse()

        # 物理量表示用テキスト
        self.info_text = OnscreenText(
            text='', pos=(-1.3, 0.9), scale=0.05,
            fg=(1, 1, 1, 1), align=TextNode.ALeft
        )

        self.sim_time = 0.0
        self.frame_count = 0
        self.taskMgr.add(self.update, "updateTask")

    def create_galaxies(self):
        # 銀河1（左側、青色）
        for _ in range(25):
            r = random.uniform(2, 15)
            theta = random.uniform(0, 2 * math.pi)
            x = -30.0 + r * math.cos(theta)
            y = r * math.sin(theta)
            z = random.uniform(-2, 2)
            orbital_speed = 1.5 * math.sqrt(1.0 / max(r, 1.0))
            vx = 2.0 - orbital_speed * math.sin(theta)
            vy = orbital_speed * math.cos(theta)
            vz = random.uniform(-0.2, 0.2)
            mass = random.uniform(0.5, 1.5)
            self.bodies.append(Body(x, y, z, vx, vy, vz, mass, (0.3, 0.5, 1.0, 1.0)))

        # 銀河2（右側、赤色）
        for _ in range(25):
            r = random.uniform(2, 15)
            theta = random.uniform(0, 2 * math.pi)
            x = 30.0 + r * math.cos(theta)
            y = r * math.sin(theta)
            z = random.uniform(-2, 2)
            orbital_speed = 1.5 * math.sqrt(1.0 / max(r, 1.0))
            vx = -2.0 - orbital_speed * math.sin(theta)
            vy = orbital_speed * math.cos(theta)
            vz = random.uniform(-0.2, 0.2)
            mass = random.uniform(0.5, 1.5)
            self.bodies.append(Body(x, y, z, vx, vy, vz, mass, (1.0, 0.3, 0.3, 1.0)))

    def create_particle_mesh(self):
        fmt = GeomVertexFormat.getV3c4()
        vdata = GeomVertexData('bodies', fmt, Geom.UHDynamic)
        vertex = GeomVertexWriter(vdata, 'vertex')
        color = GeomVertexWriter(vdata, 'color')
        for body in self.bodies:
            vertex.addData3(*body.pos)
            color.addData4(*body.color)
        points = GeomPoints(Geom.UHDynamic)
        points.addNextVertices(len(self.bodies))
        geom = Geom(vdata)
        geom.addPrimitive(points)
        node = GeomNode('bodies_node')
        node.addGeom(geom)
        return node

    def create_trail_mesh(self):
        fmt = GeomVertexFormat.getV3c4()
        vdata = GeomVertexData('trails', fmt, Geom.UHDynamic)
        vdata.setNumRows(0)
        geom = Geom(vdata)
        node = GeomNode('trails_node')
        node.addGeom(geom)
        return node

    def calculate_forces(self):
        """全天体間の重力を計算"""
        n = len(self.bodies)
        forces = np.zeros((n, 3), dtype=np.float64)
        for i in range(n):
            for j in range(i + 1, n):
                diff = self.bodies[j].pos - self.bodies[i].pos
                r_squared = np.dot(diff, diff) + self.softening ** 2
                r = np.sqrt(r_squared)
                f_mag = self.G * self.bodies[i].mass * self.bodies[j].mass / r_squared
                force = f_mag * diff / r
                forces[i] += force
                forces[j] -= force
        return forces

    def update_physics(self):
        """オイラー法による時間積分"""
        forces = self.calculate_forces()
        for i, body in enumerate(self.bodies):
            acc = forces[i] / body.mass
            body.vel += acc * self.dt
            body.pos += body.vel * self.dt

    def calculate_energy(self):
        kinetic = 0.0
        for body in self.bodies:
            kinetic += 0.5 * body.mass * np.dot(body.vel, body.vel)
        potential = 0.0
        n = len(self.bodies)
        for i in range(n):
            for j in range(i + 1, n):
                diff = self.bodies[j].pos - self.bodies[i].pos
                r = np.sqrt(np.dot(diff, diff) + self.softening ** 2)
                potential -= self.G * self.bodies[i].mass * self.bodies[j].mass / r
        return kinetic, potential, kinetic + potential

    def calculate_momentum(self):
        p = np.zeros(3)
        for body in self.bodies:
            p += body.mass * body.vel
        return p

    def update_mesh(self):
        geom = self.particle_mesh.modifyGeom(0)
        vdata = geom.modifyVertexData()
        vertex = GeomVertexWriter(vdata, 'vertex')
        for body in self.bodies:
            vertex.setData3(*body.pos)

    def update_trails(self):
        geom = self.trail_mesh.modifyGeom(0)
        vdata = geom.modifyVertexData()
        total_vertices = sum((len(b.trail) - 1) * 2 for b in self.bodies if len(b.trail) >= 2)
        vdata.setNumRows(total_vertices)
        vertex = GeomVertexWriter(vdata, 'vertex')
        color = GeomVertexWriter(vdata, 'color')
        geom.clearPrimitives()
        lines = GeomLines(Geom.UHDynamic)
        idx = 0
        for body in self.bodies:
            if len(body.trail) < 2:
                continue
            for i in range(len(body.trail) - 1):
                alpha = (i + 1) / len(body.trail) * 0.7
                vertex.setData3(*body.trail[i])
                color.setData4(body.color[0], body.color[1], body.color[2], alpha)
                vertex.setData3(*body.trail[i + 1])
                color.setData4(body.color[0], body.color[1], body.color[2], alpha)
                lines.addVertices(idx, idx + 1)
                idx += 2
        geom.addPrimitive(lines)

    def update_info_text(self):
        kinetic, potential, total = self.calculate_energy()
        p_mag = np.linalg.norm(self.calculate_momentum())
        self.info_text.setText(
            f"Time: {self.sim_time:.1f}\n"
            f"Bodies: {len(self.bodies)}\n"
            f"Kinetic E: {kinetic:.2f}\n"
            f"Potential E: {potential:.2f}\n"
            f"Total E: {total:.2f}\n"
            f"Momentum: {p_mag:.2f}"
        )

    def update(self, task):
        # 物理演算を5回実行（安定性向上）
        for _ in range(5):
            self.update_physics()
            self.sim_time += self.dt

        # 5フレームごとに軌跡を記録
        if self.frame_count % 5 == 0:
            for body in self.bodies:
                body.add_trail_point()

        self.update_mesh()
        self.update_trails()
        self.update_info_text()
        self.frame_count += 1
        return Task.cont


app = GalaxyCollision()
app.run()

群れ行動シミュレーション（Boidsアルゴリズム）

このプログラムが描くもの

立方体空間（一辺52）の中に、500個の三角錐（個体）がランダムな位置と速度で配置される。各個体は毎フレーム、自分から半径10以内の仲間を見て、次の3つのルールで進路を決める。

• 分離：近すぎる仲間からは離れる
• 整列：仲間の進む向きに合わせる
• 結合：仲間の重心に向かう

これはCraig Reynoldsが1987年に提案したBoidsアルゴリズムである。各個体は進行方向に先端を向けて描画される。空間の端から出た個体は反対側から再出現する。

画面左上には5つの指標（平均速度、速度のばらつき、平均近傍数、空間の広がり、整列スコア）が表示され、それぞれが正常範囲内かどうかも判定される。

3次元グラフィックスとして見るべき点

本プログラムは「複雑な形状を作る技術」ではなく、「たくさんの単純な物体を動かして全体として何かを表現する」技術の例である。映画やゲームで群衆や魚群を描くときの基礎になっている。各個体は単純な三角錐だが、500個が協調して動くと、1つの生き物のような印象を与える。

演習

《手順》

1. 共通の準備の手順でコードを実行する。最初はランダムな500個体が、数秒～十数秒で群れ始める
2. 群れができる前（最初の数秒）と、群れができた後を比較観察する
3. マウス操作で視点を動かし、横からと進行方向の正面からの両方で観察する
4. 画面左上の5つの指標を確認（実行直後・1分後・3分後）する。すべて「OK」と表示されていれば、群れが健全に形成されている
5. 群れが分裂したり再合流したりする様子を待ち、その瞬間に指標がどう変化するか観察する

《ヒント》

• 各個体は半径10以内しか「見えない」。空間サイズは一辺52なので、全体を見渡すことはできない。それでも群れができることが、この演習の核心である
• 「整列スコア」は0～1の値。1に近いほど全員が同じ向きに進んでいる。0.3を超えれば群れと言える
• 群れは1つにまとまるとは限らない。複数の群れに分かれることも、再び合流することもある
• 視点を引いて全体を見るときと、近づいて1個体を追うときでは、見え方が大きく変わる

《考察ポイント》

• 個と全体：1個体だけ目で追ってみる。その動きは規則的に見えるか、ランダムに見えるか。同じ動きを500個まとめて見たときの印象との違いを言葉にする
• ルールと見た目の関係：3つのルール（離れる・合わせる・近づく）のうち、もし1つだけ抜けたら見た目はどう変わると予想されるか
• 三角錐という選択：個体を球（向きが分からない）ではなく三角錐（先端で向きが分かる）にしている。これにより群れの「流れ」が見える。形の選択が、伝わる情報量を決める
• 応用：このアルゴリズムは映画の群衆シーン、戦略ゲームの大軍勢、水族館展示の魚群シミュレータなどに使われている。1個ずつ手で振り付けることが不可能なシーンを成立させている

# コード2：群れ行動シミュレーション（Boidsアルゴリズム）
from direct.showbase.ShowBase import ShowBase
from direct.task import Task
from panda3d.core import GeomVertexFormat, GeomVertexData, GeomVertexWriter
from panda3d.core import Geom, GeomTriangles, GeomNode, Mat4
from panda3d.core import LPoint3f, TextNode
from direct.gui.OnscreenText import OnscreenText
import numpy as np


class BoidsSimulation(ShowBase):
    def __init__(self):
        ShowBase.__init__(self)

        # 空間サイズと個体数
        self.space_size = 26.0
        self.num_boids = 500

        # Boidパラメータ
        self.max_speed = 8.0
        self.max_force = 0.4
        self.perception_radius = 10.0

        # 全個体の状態
        self.positions = np.random.uniform(-6.5, 6.5, (self.num_boids, 3))
        self.velocities = np.random.uniform(-3, 3, (self.num_boids, 3))

        # 各個体の可視化ノード（三角錐）
        self.boid_nodes = []
        for _ in range(self.num_boids):
            node = self.create_boid_mesh()
            self.boid_nodes.append(self.render.attachNewNode(node))

        # 指標表示用テキスト
        self.info_text = OnscreenText(
            text='', pos=(-1.3, 0.9), scale=0.05,
            fg=(1, 1, 1, 1), align=TextNode.ALeft
        )

        # カメラ設定
        self.disableMouse()
        self.camera.setPos(0, -120, 60)
        self.camera.lookAt(0, 0, 0)
        self.camLens.setFov(90)
        mat = Mat4(self.camera.getMat())
        mat.invertInPlace()
        self.mouseInterfaceNode.setMat(mat)
        self.enableMouse()

        self.frame_count = 0
        self.taskMgr.add(self.update, "updateTask")
        self.prev_time = globalClock.getFrameTime()

    def create_boid_mesh(self):
        """三角錐メッシュを作成"""
        fmt = GeomVertexFormat.getV3c4()
        vdata = GeomVertexData('boid', fmt, Geom.UHStatic)
        vertex = GeomVertexWriter(vdata, 'vertex')
        color = GeomVertexWriter(vdata, 'color')
        for v in [(0, 1, 0), (-0.3, -0.5, 0.3), (0.3, -0.5, 0.3), (0, -0.5, -0.3)]:
            vertex.addData3(*v)
            color.addData4(0.8, 0.3, 0.2, 1.0)
        tris = GeomTriangles(Geom.UHStatic)
        tris.addVertices(0, 1, 2)
        tris.addVertices(0, 2, 3)
        tris.addVertices(0, 3, 1)
        tris.addVertices(1, 3, 2)
        geom = Geom(vdata)
        geom.addPrimitive(tris)
        node = GeomNode('boid_node')
        node.addGeom(geom)
        return node

    def calculate_all_forces(self):
        """全個体の力を計算（分離・整列・結合）"""
        diff = self.positions[:, np.newaxis, :] - self.positions[np.newaxis, :, :]
        distances = np.linalg.norm(diff, axis=2)
        np.fill_diagonal(distances, np.inf)
        in_range = distances < self.perception_radius
        neighbor_counts = np.sum(in_range, axis=1)
        counts_safe = np.where(neighbor_counts > 0, neighbor_counts, 1)

        # 分離: 距離の2乗で割った反発ベクトル
        dist_safe = np.where(distances > 0, distances, 1)
        sep_w = np.where(in_range, 1.0 / (dist_safe ** 2), 0)
        separation = np.sum(diff * sep_w[:, :, np.newaxis], axis=1) / counts_safe[:, np.newaxis]

        # 整列: 近傍の速度の平均
        v_masked = self.velocities[np.newaxis, :, :] * in_range[:, :, np.newaxis]
        alignment = np.sum(v_masked, axis=1) / counts_safe[:, np.newaxis]

        # 結合: 近傍の重心への方向
        p_masked = self.positions[np.newaxis, :, :] * in_range[:, :, np.newaxis]
        cohesion = np.sum(p_masked, axis=1) / counts_safe[:, np.newaxis] - self.positions

        force = self.steer(separation) * 1.8 + self.steer(alignment) * 1.5 + self.steer(cohesion) * 1.8
        return np.where(neighbor_counts[:, np.newaxis] > 0, force, 0)

    def steer(self, desired):
        """目標方向へのステアリング力"""
        norms = np.linalg.norm(desired, axis=1, keepdims=True)
        norms_safe = np.where(norms > 0, norms, 1)
        target = (desired / norms_safe) * self.max_speed
        steering = target - self.velocities
        s_norms = np.linalg.norm(steering, axis=1, keepdims=True)
        s_norms_safe = np.where(s_norms > 0, s_norms, 1)
        steering = np.where(s_norms > self.max_force,
                            steering / s_norms_safe * self.max_force, steering)
        return np.where(norms > 0, steering, 0)

    def update_physics(self, dt):
        forces = self.calculate_all_forces()
        self.velocities += forces * dt

        # 速度を最大値に制限
        speeds = np.linalg.norm(self.velocities, axis=1, keepdims=True)
        speeds_safe = np.where(speeds > 0, speeds, 1)
        self.velocities = np.where(speeds > self.max_speed,
                                    self.velocities / speeds_safe * self.max_speed,
                                    self.velocities)

        self.positions += self.velocities * dt

        # トーラス型境界
        self.positions = np.where(self.positions > self.space_size,
                                   self.positions - 2 * self.space_size, self.positions)
        self.positions = np.where(self.positions < -self.space_size,
                                   self.positions + 2 * self.space_size, self.positions)

    def calculate_metrics(self):
        """5つの指標を計算"""
        speeds = np.linalg.norm(self.velocities, axis=1)
        avg_speed = np.mean(speeds)
        speed_variance = np.var(speeds)

        diff = self.positions[:, np.newaxis, :] - self.positions[np.newaxis, :, :]
        distances = np.linalg.norm(diff, axis=2)
        np.fill_diagonal(distances, np.inf)
        in_range = distances < self.perception_radius
        avg_neighbors = np.mean(np.sum(in_range, axis=1))

        center = np.mean(self.positions, axis=0)
        spatial_distribution = np.mean(np.linalg.norm(self.positions - center, axis=1))

        norms = np.linalg.norm(self.velocities, axis=1, keepdims=True)
        norms_safe = np.where(norms > 0, norms, 1)
        normalized = self.velocities / norms_safe
        alignment_score = np.linalg.norm(np.sum(normalized, axis=0)) / self.num_boids

        return {
            'avg_speed': avg_speed,
            'speed_variance': speed_variance,
            'avg_neighbors': avg_neighbors,
            'spatial_distribution': spatial_distribution,
            'alignment_score': alignment_score
        }

    def evaluate_quality(self, m):
        """各指標が正常範囲内かを判定"""
        checks = []
        checks.append("Speed: OK" if 2.0 < m['avg_speed'] < 8.0
                      else f"Speed: ABNORMAL ({m['avg_speed']:.2f})")
        checks.append("Speed variance: OK" if m['speed_variance'] < 3.0
                      else f"Speed variance: HIGH ({m['speed_variance']:.2f})")
        checks.append("Neighbors: OK" if 3.0 < m['avg_neighbors'] < 30.0
                      else f"Neighbors: ABNORMAL ({m['avg_neighbors']:.2f})")
        checks.append("Distribution: OK" if 3.0 < m['spatial_distribution'] < 20.0
                      else f"Distribution: ABNORMAL ({m['spatial_distribution']:.2f})")
        checks.append("Alignment: OK" if m['alignment_score'] > 0.3
                      else f"Alignment: LOW ({m['alignment_score']:.2f})")
        return checks

    def update_visuals(self):
        """各個体を位置に移動し、進行方向に向ける"""
        speeds = np.linalg.norm(self.velocities, axis=1)
        for i in range(self.num_boids):
            self.boid_nodes[i].setPos(float(self.positions[i, 0]),
                                       float(self.positions[i, 1]),
                                       float(self.positions[i, 2]))
            if speeds[i] > 0:
                target = self.positions[i] + self.velocities[i] / speeds[i]
                self.boid_nodes[i].lookAt(LPoint3f(float(target[0]),
                                                    float(target[1]),
                                                    float(target[2])))

    def update(self, task):
        current_time = globalClock.getFrameTime()
        dt = current_time - self.prev_time
        self.prev_time = current_time

        self.update_physics(dt)
        self.update_visuals()
        self.frame_count += 1

        # 30フレームごとに指標を更新
        if self.frame_count % 30 == 0:
            m = self.calculate_metrics()
            lines = [
                f"Frame: {self.frame_count}",
                f"Boids: {self.num_boids}",
                f"Avg speed: {m['avg_speed']:.2f}",
                f"Speed variance: {m['speed_variance']:.2f}",
                f"Avg neighbors: {m['avg_neighbors']:.1f}",
                f"Spatial dist: {m['spatial_distribution']:.2f}",
                f"Alignment: {m['alignment_score']:.2f}",
                "",
                "Quality Check:"
            ] + self.evaluate_quality(m)
            self.info_text.setText('\n'.join(lines))

        return Task.cont


app = BoidsSimulation()
app.run()

流体シミュレーション（煙の拡散と上昇気流）

このプログラムが描くもの

3次元空間を63×63×88の格子に区切り、各格子点に「煙の濃さ（密度）」と「流れの速さ（速度）」を持たせる。毎フレーム、次の物理過程をこの格子上で順番に計算する。

• 浮力：濃い煙ほど上向きの力を受ける
• 拡散：濃さや速さが周囲ににじみ広がる
• 移流：流れに乗って煙が運ばれる
• 減衰：時間とともに少しずつ薄くなる

そして格子の状態とは別に、最大10000個のパーティクル（白い小さな点）が空間に漂っている。各パーティクルは、自分のいる格子の流れに乗って動き、自分のいる格子の煙が濃いほど明るく表示される。これにより、目に見えない密度場が、無数のパーティクルとして目に見える煙になる。

3次元グラフィックスとして見るべき点

銀河（点）、群れ（三角錐）に続き、煙は「形のないものをどう描くか」という別種の課題を扱う。映画やゲームでの炎・煙・霧・雲などの表現は、ほぼすべてこのパーティクルの考え方を基礎にしている。物体の形を直接モデリングするのではなく、多数の小さな要素の集まりとして見せるという発想を体験する。

演習

《手順》

1. 共通の準備の手順でコードを実行する。初回起動時はNumbaのコンパイルで数秒かかる
2. 底面中央から白い煙が立ち上り、上昇しながら横に広がる様子を、少なくとも30秒間そのまま観察する
3. マウス操作で視点を動かし、煙の柱を真横からと真上からの両方で観察する
4. 煙の立ち上り口（底面付近）と、上の方の煙を比べる。濃さ・広がり・動きの違いを確認する
5. 煙の柱が時間とともに揺らぐ・蛇行する瞬間を捕まえる

《ヒント》

• パーティクルは1つでは煙に見えない。集団として見たときに初めて煙に見える。視点を引いて全体を見るのと、近づいて1つ1つの粒を追うのとでは、見え方がまったく異なる
• 真上から見ると、底面の円形の煙源と、上空での横方向の広がりが分かりやすい
• 煙の揺らぎは、煙源にランダムな横向きの速度が加えられていることから生まれる。完全に対称な煙源では、揺らぎがなければ単調な柱になってしまう
• パーティクル数（最大10000）を見たとき、これより少なかったら煙に見えるか、多かったら何が変わるかを想像してみる

《考察ポイント》

• 面ではなく点の集まり：本物の煙には明確な「表面」がない。だから三角形メッシュ（普通の3Dモデリング）では表現しづらい。パーティクルという選択が、いかにこの題材に合っているかを言葉にしてみる
• 濃さと明るさ：パーティクル1個の色は、それがいる場所の密度に応じて明るくなる。これにより薄い煙は淡く、濃い煙は明るく見える。1個の表示ルールと全体の見え方の関係を考える
• 3つの題材の共通点：銀河（点）、群れ（三角錐）、煙（パーティクル）。3つとも「たくさんの単純な要素が、ルールに従って動くと、複雑で意味のある見え方が生まれる」という構造を持つ。これが3次元コンピュータグラフィックスの主要な表現戦略の1つである
• 応用：パーティクル系の表現は、映画の爆発・煙・魔法のエフェクト、ゲームの炎や水しぶき、気象シミュレーションの可視化などに使われている。本プログラムはその原型である

# コード3：流体シミュレーション（煙の拡散と上昇気流）
from direct.showbase.ShowBase import ShowBase
from direct.task import Task
from panda3d.core import GeomVertexFormat, GeomVertexData, GeomVertexWriter
from panda3d.core import Mat4
from panda3d.core import Geom, GeomNode, GeomPoints, RenderModeAttrib
import numpy as np
from numba import jit, prange
import random
import math


@jit(nopython=True, parallel=True, cache=True)
def advect(field, vx, vy, vz, dt, spacing, gx, gy, gz):
    """移流：速度場に沿ってフィールド値を移動（セミ・ラグランジュ法）"""
    new_field = np.zeros_like(field)
    for i in prange(1, gx - 1):
        for j in range(1, gy - 1):
            for k in range(1, gz - 1):
                # バックトレース
                x = i - dt * vx[i, j, k] / spacing * 10
                y = j - dt * vy[i, j, k] / spacing * 10
                z = k - dt * vz[i, j, k] / spacing * 10
                x = max(0.5, min(gx - 1.5, x))
                y = max(0.5, min(gy - 1.5, y))
                z = max(0.5, min(gz - 1.5, z))
                i0 = int(x); j0 = int(y); k0 = int(z)
                i1 = min(i0 + 1, gx - 1)
                j1 = min(j0 + 1, gy - 1)
                k1 = min(k0 + 1, gz - 1)
                sx = x - i0; sy = y - j0; sz = z - k0
                # 3次元線形補間
                new_field[i, j, k] = (
                    field[i0, j0, k0] * (1-sx) * (1-sy) * (1-sz) +
                    field[i1, j0, k0] * sx * (1-sy) * (1-sz) +
                    field[i0, j1, k0] * (1-sx) * sy * (1-sz) +
                    field[i0, j0, k1] * (1-sx) * (1-sy) * sz +
                    field[i1, j1, k0] * sx * sy * (1-sz) +
                    field[i1, j0, k1] * sx * (1-sy) * sz +
                    field[i0, j1, k1] * (1-sx) * sy * sz +
                    field[i1, j1, k1] * sx * sy * sz
                )
    return new_field


@jit(nopython=True, parallel=True, cache=True)
def diffuse(field, diff_rate, dt, gx, gy, gz):
    """拡散：ガウス・ザイデル法による4回反復"""
    avg_size = (gx + gy + gz) / 3
    a = dt * diff_rate * (avg_size - 2) ** 2
    new_field = field.copy()
    for _ in range(4):
        for i in prange(1, gx - 1):
            for j in range(1, gy - 1):
                for k in range(1, gz - 1):
                    new_field[i, j, k] = (
                        field[i, j, k] + a * (
                            new_field[i-1, j, k] + new_field[i+1, j, k] +
                            new_field[i, j-1, k] + new_field[i, j+1, k] +
                            new_field[i, j, k-1] + new_field[i, j, k+1]
                        )
                    ) / (1 + 6 * a)
    return new_field


class FluidSimulation(ShowBase):
    def __init__(self):
        ShowBase.__init__(self)
        self.setBackgroundColor(0, 0, 0)

        # 格子サイズ
        self.gx = 63
        self.gy = 63
        self.gz = 88
        self.spacing = 0.5

        # 密度場と速度場
        shape = (self.gx, self.gy, self.gz)
        self.density = np.zeros(shape, dtype=np.float64)
        self.velocity_x = np.zeros(shape, dtype=np.float64)
        self.velocity_y = np.zeros(shape, dtype=np.float64)
        self.velocity_z = np.zeros(shape, dtype=np.float64)

        # シミュレーションパラメータ
        self.diffusion = 0.0002
        self.viscosity = 0.0001
        self.buoyancy = 0.3
        self.dissipation = 0.995

        # パーティクル
        self.particles = []
        self.max_particles = 10000

        # パーティクルメッシュ
        self.particle_mesh = self.create_particle_mesh()
        node = self.render.attachNewNode(self.particle_mesh)
        node.setPos(-self.gx * self.spacing / 2, -self.gy * self.spacing / 2, 0)
        node.setRenderMode(RenderModeAttrib.MPoint, 8.0)

        # カメラ設定
        self.disableMouse()
        self.camera.setPos(0, -65, 25)
        self.camera.lookAt(0, 0, 18)
        self.camLens.setFov(80)
        mat = Mat4(self.camera.getMat())
        mat.invertInPlace()
        self.mouseInterfaceNode.setMat(mat)
        self.enableMouse()

        self.taskMgr.add(self.update, "updateTask")
        self.prev_time = globalClock.getFrameTime()

    def create_particle_mesh(self):
        fmt = GeomVertexFormat.getV3c4()
        vdata = GeomVertexData('particles', fmt, Geom.UHDynamic)
        vertex = GeomVertexWriter(vdata, 'vertex')
        color = GeomVertexWriter(vdata, 'color')
        for _ in range(self.max_particles):
            vertex.addData3(0, 0, 0)
            color.addData4(1.0, 1.0, 1.0, 0.0)
        points = GeomPoints(Geom.UHDynamic)
        points.addNextVertices(self.max_particles)
        geom = Geom(vdata)
        geom.addPrimitive(points)
        node = GeomNode('particle_node')
        node.addGeom(geom)
        return node

    def add_smoke_source(self):
        """底面中央に煙源を追加"""
        cx = self.gx // 2
        cy = self.gy // 2
        radius = 6
        for i in range(cx - radius, cx + radius):
            for j in range(cy - radius, cy + radius):
                if math.sqrt((i - cx) ** 2 + (j - cy) ** 2) < radius:
                    for k in range(0, 3):
                        self.density[i, j, k] = min(1.0, self.density[i, j, k] + 0.3)
                        self.velocity_z[i, j, k] = 3.0
                        self.velocity_x[i, j, k] += random.uniform(-0.5, 0.5)
                        self.velocity_y[i, j, k] += random.uniform(-0.5, 0.5)

    def update_fluid(self, dt):
        """流体シミュレーションの更新"""
        dt = min(dt, 0.05)
        gx, gy, gz = self.gx, self.gy, self.gz

        # 1. 浮力
        self.velocity_z += self.buoyancy * self.density * dt

        # 2. 速度場の拡散
        self.velocity_x = diffuse(self.velocity_x, self.viscosity, dt, gx, gy, gz)
        self.velocity_y = diffuse(self.velocity_y, self.viscosity, dt, gx, gy, gz)
        self.velocity_z = diffuse(self.velocity_z, self.viscosity, dt, gx, gy, gz)

        # 3. 速度場の移流
        self.velocity_x = advect(self.velocity_x, self.velocity_x, self.velocity_y, self.velocity_z, dt, self.spacing, gx, gy, gz)
        self.velocity_y = advect(self.velocity_y, self.velocity_x, self.velocity_y, self.velocity_z, dt, self.spacing, gx, gy, gz)
        self.velocity_z = advect(self.velocity_z, self.velocity_x, self.velocity_y, self.velocity_z, dt, self.spacing, gx, gy, gz)

        # 4. 密度場の拡散と移流
        self.density = diffuse(self.density, self.diffusion, dt, gx, gy, gz)
        self.density = advect(self.density, self.velocity_x, self.velocity_y, self.velocity_z, dt, self.spacing, gx, gy, gz)

        # 5. 減衰
        self.density *= self.dissipation

        # 6. 煙源追加
        self.add_smoke_source()

    def update_particles(self):
        """パーティクルの位置と不透明度を更新"""
        gx, gy, gz = self.gx, self.gy, self.gz

        # 既存パーティクルの更新
        for particle in self.particles[:]:
            i = int(particle['x'] / self.spacing)
            j = int(particle['y'] / self.spacing)
            k = int(particle['z'] / self.spacing)
            if 0 <= i < gx - 1 and 0 <= j < gy - 1 and 0 <= k < gz - 1:
                particle['x'] += self.velocity_x[i, j, k] * 0.05
                particle['y'] += self.velocity_y[i, j, k] * 0.05
                particle['z'] += self.velocity_z[i, j, k] * 0.05
                particle['alpha'] = min(1.0, self.density[i, j, k] * 1.5)
                particle['life'] -= 1
            else:
                particle['life'] = 0
            if particle['life'] <= 0:
                self.particles.remove(particle)

        # 新しいパーティクルの追加
        if len(self.particles) < self.max_particles:
            cx = self.gx // 2
            cy = self.gy // 2
            for _ in range(min(100, self.max_particles - len(self.particles))):
                x = (cx + random.uniform(-5, 5)) * self.spacing
                y = (cy + random.uniform(-5, 5)) * self.spacing
                z = random.uniform(0, 1) * self.spacing
                self.particles.append({'x': x, 'y': y, 'z': z, 'alpha': 0.8, 'life': 450})

    def update_mesh(self):
        """パーティクルメッシュを更新"""
        geom = self.particle_mesh.modifyGeom(0)
        vdata = geom.modifyVertexData()
        vertex = GeomVertexWriter(vdata, 'vertex')
        color = GeomVertexWriter(vdata, 'color')
        for particle in self.particles:
            vertex.setData3(particle['x'], particle['y'], particle['z'])
            alpha = particle['alpha']
            brightness = 0.7 + alpha * 0.3
            color.setData4(brightness, brightness, brightness, alpha)
        # 未使用スロットは透明に
        for _ in range(len(self.particles), self.max_particles):
            vertex.setData3(0, 0, 0)
            color.setData4(0, 0, 0, 0)

    def update(self, task):
        current_time = globalClock.getFrameTime()
        dt = current_time - self.prev_time
        self.prev_time = current_time

        self.update_fluid(dt)
        self.update_particles()
        self.update_mesh()
        return Task.cont


app = FluidSimulation()
app.run()