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

前準備：環境構築とプログラムの実行方法

Python 3.12 のインストール

Pythonのインストールを行い、Pythonのプログラムを実行する環境を整える。扱う環境は、Windows搭載パソコンである。金子研究室では、Python 3.12.10を推奨する。

REM Python 3.12 をシステム領域にインストール
winget install --id Python.Python.3.12 -e --scope machine --silent --accept-source-agreements --accept-package-agreements --override "/quiet InstallAllUsers=1 PrependPath=1 Include_test=0 Include_pip=1 Include_launcher=1 InstallLauncherAllUsers=1 TargetDir=\"C:\Program Files\Python312\""

REM Python と Scripts を PATH 先頭に追加
powershell -NoProfile -Command "$p='C:\Program Files\Python312'; $s=\"$p\Scripts\"; $c=[Environment]::GetEnvironmentVariable('Path','Machine'); if((Test-Path $p) -and (';'+$c+';' -notlike \"*;$p;*\") -and (';'+$c+';' -notlike \"*;$s;*\")){[Environment]::SetEnvironmentVariable('Path',\"$p;$s;$c\",'Machine')}"

python --version

Python の開発環境 Visual Studio Code のインストールと Python 用の設定

Python の開発環境Visual Studio Code（プログラムを編集するソフトウェア。以下、VS Code）を整える。

REM ============================================================
REM Microsoft Visual Studio Code
REM ============================================================
winget uninstall -e --id Microsoft.VisualStudioCode --silent --disable-interactivity --accept-source-agreements
rmdir /s /q C:\ProgramData\vscode-extensions 2>nul
rmdir /s /q "%APPDATA%\Code" 2>nul
rmdir /s /q "%USERPROFILE%\.vscode" 2>nul
rmdir /s /q "%LOCALAPPDATA%\Microsoft\vscode-update" 2>nul

REM VS Code をシステム領域に新規インストール
winget install --scope machine --id Microsoft.VisualStudioCode -e --silent --accept-source-agreements --accept-package-agreements

REM 全ユーザー共有の拡張機能フォルダ
mkdir C:\ProgramData\vscode-extensions 2>nul
icacls "C:\ProgramData\vscode-extensions" /grant "Everyone:(OI)(CI)M" /T

REM スタートメニューのショートカットを --extensions-dir 付きで再作成
rmdir /s /q "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code" 2>nul
del "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk" 2>nul
powershell -NoProfile -Command "$s=New-Object -ComObject WScript.Shell; $lnk=$s.CreateShortcut('C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk'); $lnk.TargetPath='C:\Program Files\Microsoft VS Code\Code.exe'; $lnk.Arguments='--extensions-dir \"C:\ProgramData\vscode-extensions\"'; $lnk.Save()"
REM ショートカットの検証
powershell -NoProfile -Command "$s=New-Object -ComObject WScript.Shell; $lnk=$s.CreateShortcut('C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk'); Write-Host 'TargetPath:' $lnk.TargetPath; Write-Host 'Arguments:' $lnk.Arguments"

REM ファイル / フォルダ右クリックの「Code で開く」を登録
reg add "HKLM\SOFTWARE\Classes\*\shell\VSCode\command" /ve /d "\"C:\Program Files\Microsoft VS Code\Code.exe\" --extensions-dir \"C:\ProgramData\vscode-extensions\" \"%1\"" /f
reg add "HKLM\SOFTWARE\Classes\Directory\shell\VSCode\command" /ve /d "\"C:\Program Files\Microsoft VS Code\Code.exe\" --extensions-dir \"C:\ProgramData\vscode-extensions\" \"%1\"" /f
reg add "HKLM\SOFTWARE\Classes\Directory\Background\shell\VSCode\command" /ve /d "\"C:\Program Files\Microsoft VS Code\Code.exe\" --extensions-dir \"C:\ProgramData\vscode-extensions\" \"%V\"" /f

REM --extensions-dir 付きで起動する code.cmd ラッパを作成
REM （%* を echo で書くと対話的 cmd で失われるため、PowerShell で [char]37+'*' を書き出す）
powershell -NoProfile -Command "$pct=[char]37; $q=[char]34; $c='@echo off'+[char]13+[char]10+$q+'C:\Program Files\Microsoft VS Code\bin\code.cmd'+$q+' --extensions-dir '+$q+'C:\ProgramData\vscode-extensions'+$q+' '+$pct+'*'+[char]13+[char]10; [IO.File]::WriteAllText('C:\ProgramData\vscode-extensions\vscode.cmd',$c,[Text.Encoding]::ASCII)"

REM 拡張機能のインストール
set "CODE=C:\Program Files\Microsoft VS Code\bin\code.cmd"
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --uninstall-extension GitHub.copilot
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --uninstall-extension GitHub.copilot-chat
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension ms-python.python
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension ms-python.vscode-pylance
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension ms-python.debugpy
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension MS-CEINTL.vscode-language-pack-ja
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension saoudrizwan.claude-dev
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension rust-lang.rust-analyzer
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension tamasfe.even-better-toml
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension anthropic.claude-code
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension almenon.arepl
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --list-extensions --show-versions
echo === セットアップ完了 ===

Python プログラム実行手順

必要なライブラリのインストール

3つの演習はすべてPanda3Dを使う。Panda3Dは3次元コンピュータグラフィックスのためのライブラリである。インストールは1回だけ行えばよい。演習3はnumbaも使う。numbaはCPU上でPythonの数値計算を高速化するライブラリで、GPUは不要である。GPU搭載機でもCPUのみの機でも同じ手順で動作する。

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

pip install --no-user -U panda3d numba

プログラムの実行

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

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

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

python a.py

共通のマウス操作

3つの演習は同じマウス操作で視点を動かせる。Panda3Dの既定のカメラ操作（トラックボール操作）は次のとおりである。

• 左ドラッグ：視点を左右・上下に平行移動
• 右ドラッグ：ズームイン・ズームアウト（前後移動）
• 中ドラッグ：原点を中心に視点を回転

3次元コンピュータグラフィックスでは、見る側が視点を変えて対象を観察できる。各演習では視点を変えながら観察する。

プログラムの終了

表示されたウィンドウの閉じるボタンを押すか、ウィンドウを選択した状態でEscキーを押す。

共通の前提：近似計算と見え方

観察の前に、3つの演習に共通する次の前提を理解しておく。これを知らないと、正常な挙動を不具合と誤解しやすい。

(1) すべて近似計算である。3つの演習は、連続的な物理現象を細かい時間刻みと格子に区切って近似計算している（時間積分や反復計算など）。厳密解ではないため、計算を続けると、本来一定に保たれるはずの物理量が少しずつずれていく。これは数値シミュレーションの性質であり、操作ミスや不具合ではない。

(2) 表示される物理量は「急変しないか」を見る。画面に表示されるエネルギーや運動量は、理想的には一定に保たれる量である。ただし(1)の理由で緩やかにずれる。観察では、値が完全に一定かではなく、短時間で桁が変わるような急変をしていないかを確認する。急変は計算が破綻した兆候だが、緩やかなずれは正常である。

(3) 濃淡や半透明は描画上の表現である。各演習では、色の濃さ・薄さ（半透明）で軌跡の新旧や煙の濃さを表す。点や粒は奥行きに関わらず画面上で同じ大きさで描かれる（遠くの点が近くの点より小さく見えることはない）。

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

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

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

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

形状のモデルを用意せず、点と線だけで天文現象を可視化している。

背景となる考え方

万有引力の法則により、2つの物体間の引力は両者の質量の積に比例し、距離の2乗に反比例する。N個の天体があれば、すべてのペアについてこの計算を行い、各天体に働く力の合計から加速度を求め、速度と位置を時間刻みだけ進める。これを毎フレーム繰り返すと運動が時間発展する。本プログラムはこの時間積分にオイラー法（最も基本的な近似計算）を用いる。

外力が働かない孤立系では、全エネルギーと全運動量が理想的には保存される。ただしオイラー法は近似のため、これらの量は緩やかにずれていく（共通の前提(1)(2)）。画面左上の表示で確認すべきは、値が完全に一定かではなく、短時間で急変せず緩やかな変化にとどまっているかである。

演習1：銀河衝突の観察

手順

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

ヒント

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

考察ポイント

• 軌跡の半透明表示：古い軌跡ほど薄くなる表示から、運動の方向と速さが読み取れる。一定の濃さだった場合との見え方の違いを考える。
• 3次元で見る効果：2次元の図と、視点を変えられる3次元表示とで、衝突の理解にどんな違いがあるかを考える。
• 物理量の変化：Total EとMomentumが急変せず緩やかにとどまっていたかを確認する。短時間で大きく変動した場合は、天体どうしが極端に接近して計算が破綻しかけたことが原因である。
• ルールから生まれる形：衝突後の銀河の形を指定するコードは存在しない。引力の法則だけで複雑な形状が生まれることを確認する。

# コード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 panda3d.core import TransparencyAttrib
from direct.gui.OnscreenText import OnscreenText
import numpy as np
import random
import math
import sys


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.accept("escape", sys.exit)

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

        # 天体リスト
        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()
        trail_node = self.render.attachNewNode(self.trail_mesh)
        trail_node.setTransparency(TransparencyAttrib.MAlpha)
        trail_node.setDepthWrite(False)

        # カメラ設定
        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アルゴリズム）

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

中心から各軸方向に±26（全長52）の立方体空間に、500個の三角錐（個体）を配置する。個体は最初から空間全体に散らばっているのではなく、中心付近の狭い範囲（中心から±6.5、全長13）に集めて初期配置する。各個体は毎フレーム、自分から半径10以内の仲間を見て、次の3つのルールで進路を決める。半径10は初期配置の広がり（全長13）と同程度なので、開始直後から多くの個体が互いを認識でき、すぐに群れ始める。

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

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

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

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

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

演習2：群れの形成の観察

手順

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

ヒント

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

考察ポイント

• 個と全体：1個体だけを目で追い、その動きが規則的かランダムかを見る。500個まとめて見たときの印象との違いを言葉にする。
• ルールと見た目：3つのルール（離れる・合わせる・近づく）のうち1つが抜けたら見た目がどう変わるかを予想する。
• 三角錐の効果：球（向きが分からない）ではなく三角錐（先端で向きが分かる）にすることで、群れの流れが見えることを確認する。
• 指標の変化：群れの分裂・合流の瞬間に、整列スコアや空間の広がりがどう変化するかを読み取る。

# コード2：群れ行動シミュレーション（Boidsアルゴリズム）
from direct.showbase.ShowBase import ShowBase
from direct.showbase.ShowBaseGlobal import globalClock
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
import sys


class BoidsSimulation(ShowBase):
    def __init__(self):
        ShowBase.__init__(self)
        self.accept("escape", sys.exit)

        # 空間サイズ（中心から各軸±space_size、全長は2倍）と個体数
        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

        # 全個体の状態（中心付近の狭い範囲±6.5に初期配置）
        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の格子に区切り、各格子点に煙の濃さ（密度）と流れの速さ（速度）を持たせる。毎フレーム、次の物理過程を格子上で順番に計算する。計算はいずれも近似であり、初回起動時はNumbaによる機械語への変換（コンパイル）に数秒かかる。

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

格子の状態とは別に、最大10000個のパーティクル（白い小さな点）が空間に漂っている。各パーティクルは、自分のいる格子の流れに乗って動き、自分のいる格子の煙が濃いほど明るく、かつ不透明に表示される。この濃淡は半透明合成で描かれる（共通の前提(3)）。これにより、目に見えない密度場が、無数のパーティクルとして目に見える煙になる。

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

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

演習3：煙の挙動の観察

手順

1. 共通の準備の手順でコードを実行する。初回起動時はNumbaのコンパイルで数秒かかる（2回目以降はキャッシュにより短縮される）。
2. 底面中央から白い煙が立ち上り、上昇しながら横に広がる様子を、少なくとも30秒間観察する。
3. マウスの中ドラッグで視点を回転し、煙の柱を真横からと真上からの両方で観察する。
4. 煙の立ち上り口（底面付近）と上の方の煙を比べ、濃さ・広がり・動きの違いを確認する。
5. 煙の柱が揺らぐ・蛇行する瞬間を捕まえる。
6. 観察を終えたら、ウィンドウの閉じるボタンを押すかEscキーで終了する。

ヒント

• パーティクルは1つでは煙に見えず、集団として見たときに煙に見える。視点を引いて全体を見るのと、近づいて1つ1つの粒を追うのとで見え方が変わる。
• 真上から見ると、底面の円形の煙源と上空での横方向の広がりが分かる。
• 煙の揺らぎは、煙源に毎フレーム小さなランダムの横向き速度が加えられていることから生まれる。横向き速度がなく煙源が完全に対称だったら、煙はまっすぐな柱になる。この非対称性が揺らぎと蛇行を生む。

考察ポイント

• 面ではなく点の集まり：本物の煙には明確な表面がなく、三角形メッシュでは表現しづらい。パーティクルがこの題材に合っている理由を言葉にする。
• 濃さと明るさ：パーティクル1個の色と不透明度は、その場所の密度で決まる。薄い煙は淡く、濃い煙は明るく見える。1個の表示ルールと全体の見え方の関係を考える。
• 3つの題材の共通点：銀河（点）、群れ（三角錐）、煙（パーティクル）はいずれも、多数の単純な要素がルールに従って動くと複雑な見え方が生まれるという構造を持つ。
• 立ち上り口と上方の違い：底面付近と上方とで煙の濃さ・広がりがどう違うかを読み取り、拡散と減衰の効果を確認する。

# コード3：流体シミュレーション（煙の拡散と上昇気流）
from direct.showbase.ShowBase import ShowBase
from direct.showbase.ShowBaseGlobal import globalClock
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
from panda3d.core import TransparencyAttrib
import numpy as np
from numba import jit, prange
import random
import math
import sys


@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, 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 range(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.accept("escape", sys.exit)

        # 格子サイズ
        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)
        node.setTransparency(TransparencyAttrib.MAlpha)
        node.setDepthWrite(False)

        # カメラ設定
        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. 速度場の移流（移流前の速度場を退避し、3成分を同じ速度場で移流する）
        vx_old = self.velocity_x.copy()
        vy_old = self.velocity_y.copy()
        vz_old = self.velocity_z.copy()
        self.velocity_x = advect(self.velocity_x, vx_old, vy_old, vz_old, dt, self.spacing, gx, gy, gz)
        self.velocity_y = advect(self.velocity_y, vx_old, vy_old, vz_old, dt, self.spacing, gx, gy, gz)
        self.velocity_z = advect(self.velocity_z, vx_old, vy_old, vz_old, 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()