7. バネ・質点系の力学シミュレーション

前章までの運動学シミュレーションでは，質点どうしの間には相互作用がありませんでした．この章では，力学シミュレーションを題材として，オブジェクト間にもう少し複雑な相互作用があるプログラムの実際を学びます．

難しそうに聞こえるかもしれませんが，これまで学んできたクラスの継承やオブジェクトの合成などの考え方を活用することで，単純な要素の組み合わせで複雑なシステムを構築できることを見ていきます．

  • 左クリック: 質点を生成
  • 右クリック: 固定質点を生成
  • スペースキーを押しながらクリック: バネでつながれた質点を生成
  • ESCキー: 終了

7.1. まず試してみる

この章で扱う力学シミュレーションは，質点，バネ，質点間の衝突，質点と壁の衝突を扱います．これらは spring_mass モジュールとして利用できるようにします．

spring_mass モジュールを作る話は後回しにして，まずはユーザとして使ってみるところから始めます．

これまで同様 projects フォルダの中に spring_mass というフォルダを作成し，VSCode でそのフォルダを開いた状態にしてください．そこに spring_mass.py というファイルを作成して，この章の最後の spring_mass.py (ver 18.0) のコードをコピーペーストして保存しておいてください．この時点で VSCode のタイトルバーには spring_mass.py - spring_mass - VSCodium (Computer Seminar I) と表示されているはずです．

次に，spring_mass モジュールを使う側のプログラムとして simple_main.py ファイルを作成し，以下の内容を入力してください．

simple_main.py (ver 1.0)
 1import pygame
 2import spring_mass as spm
 3
 4
 5def main():
 6    pygame.init()
 7    width, height = 600, 400
 8    screen = pygame.display.set_mode((width, height))
 9    clock = pygame.time.Clock()
10
11    world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.5))
12    actor_list = []
13
14
15    while True:
16        frames_per_second = 60
17        clock.tick(frames_per_second)
18
19        should_quit = False
20        for event in pygame.event.get():
21            if event.type == pygame.QUIT:
22                should_quit = True
23            elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
24                should_quit = True
25        if should_quit:
26            break
27
28        for a in actor_list:
29            a.update()
30        actor_list[:] = [a for a in actor_list if a.is_alive]
31
32        screen.fill(pygame.Color("black"))
33        for a in actor_list:
34            a.draw(screen)
35        pygame.display.update()
36
37    pygame.quit()
38
39
40if __name__ == "__main__":
41    main()

基本構造は particle.py の ver 1.0 と一緒です．

2行目で spring_mass モジュールを import しています．名前が長いので spm という名前で使えるようにしています．

以降は main 関数だけのプログラムで，いつも通りの空っぽの雛形に， World 型オブジェクトと actor_list というリストを用意しています (11-12行目)．

World の役割はこれまでと同じですが，初期化の際の引数の仕様が少し変わっています．第1引数には画面の幅と高さのタプルを渡します．第2引数に単位時間 dt，第3引数は重力加速度ベクトルをタプルで渡します (これまでは y 成分だけのスカラ値でした)．

この章では，質点やバネなどの構成要素を総称して Actor と呼ぶことにします．これまでの particle_list に相当するものを一般化して actor_list という名前にしています．

28-30行目で actor_list 内の Actor が update され，is_alive でないものは削除されます．33-34行目で各 Actor が draw されます．

真っ黒な画面が出て，Esc キーで終了できれば OK です．

Source Control サイドバーで git リポジトリを初期化し，simple_main.py をコミットしておいてください． spring_mass.py の方も一緒にコミットしてもよいですし，しなくてもよいです．

7.1.1. 質点を配置

質点を1個生成します．前章までの Particle は属性として位置，速度，半径 (と色) を持っていましたが，この章では力学的な相互作用を実現するため，質量，粘性抵抗係数，反発係数を持たせます．これらは初期化時に引数として渡すことで設定します．力学特性を持つという点で前章までとは大きく変わるので，クラス名も新たに PointMass としました．

World 型オブジェクトへの参照を持つのもこれまでと同じです．一つ大きく変わるのは，描画処理を Strategy パターンで入れ替え可能になっている点です．spring_mass モジュールに CircleDrawer という円を描画するオブジェクトのクラスが定義されています．色と線幅 (0 なら塗りつぶし) を指定して生成します．

描画処理を入れ替え可能とすることで，単に円を描くだけでなく，もっと凝ったグラフィクスを描画したり，画像を貼り付けたりすることができるようになります．具体例は後で見ることにします．

simple_main.py (ver 2.0)
11    world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.5))
12    actor_list = []
13
14    pos, vel = (100, 200), (10, -15)
15    radius, mass = 10, 10
16    viscous_damping = 0.01
17    restitution = 0.95
18    circle_drawer = spm.CircleDrawer("green", 0)
19    actor_list.append(spm.PointMass(pos, vel, world, radius, mass,
20                                    viscous_damping, restitution, circle_drawer))
21
22    while True:

18行目では緑色，塗りつぶしの circle_drawer を生成し，19行目の PointMass 型オブジェクトの生成時に渡しています．生成された PointMass 型オブジェクトは actor_list に append されています．

なお，radius 以降の引数にはデフォルト値が設定されています．具体的な値については spring_mass.py 内の定義を見てください．

7.1.2. 壁・床を配置

particle.py では壁で跳ね返るのは Particle 自身の責務 (正確には Particle が属性として保持しているストラテジオブジェクトの責務) でしたが，この章では，質点に力を及ぼす壁 (Boundary 型のオブジェクト) を生成し，質点は力を受けて結果的に跳ね返ることにします．

Boundary 型オブジェクトを生成する際は，第1引数に法線ベクトル，第2引数に壁面上に含まれる任意の1点の座標，第3引数に World 型オブジェクト，第4 引数に衝突対象になる質点を含むリストを渡します (ここに質点以外が含まれる場合は無視されます)．法線ベクトルは，壁の中に進む向きとします．

simple_main.py (ver 3.0)
21
22    actor_list.append(spm.Boundary((1, 0), (width, 0), world, actor_list))
23    actor_list.append(spm.Boundary((0, 1), (0, height), world, actor_list))
24
25    while True:

22行目で作られるのは (width, 0) を通り右向き法線を持つ境界なので，画面右端の壁になります．23行目で床を作ります．左側の壁や天井なども同様に作れます．

これらを actor_list に入れておけば，対応する境界での跳ね返りが発生します．actor_list から外せば跳ね返らず外に飛び去っていきます．いろいろ試してみてください．

ウィンドウの端じゃないところにも壁は作れますか?
作れますが「壁の中」の領域に質点があるときの運動は不自然になります．

これを解消するには，「壁の中」に一定距離以上入り込んだときは力を発生しないなどの対処が考えられます．ただしうまく調整しないと質点の速度によっては「すり抜け」が発生するかも知れません．

また，ここで作った壁には「表」と「裏」があることに注意してください．表からも裏からも衝突できる壁が欲しい場合は，表向きの壁と裏向きの壁 (それぞれ法線方向が互いと逆) の 2 枚を同じところに置くなどの対処が必要です．

7.1.3. バネを配置

バネはクラス Spring で表します．初期化時に，両端の質点 2 つと World 型オブジェクトのほか，バネ定数，自然長，描画ストラテジオブジェクトを渡します．

直線を描画するオブジェクトのクラスとして LineDrawer が用意されています．

simple_main.py (ver 4.0)
24
25    p1 = spm.PointMass((300, 100), (0, 0), world, drawer=circle_drawer)
26    p2 = spm.PointMass((400, 120), (0, 0), world, drawer=circle_drawer)
27    spring_const = 0.01
28    natural_len = 20
29    line_drawer = spm.LineDrawer("white", 3)
30    sp12 = spm.Spring(p1, p2, world, spring_const, natural_len, line_drawer)
31    actor_list += [p1, p2, sp12]
32
33    while True:

ここでは適当な 2 つの質点を生成し，バネでつないでいます．重力があるとバネ自体の動作がわかりにくいかも知れません．11行目の重力加速度を適宜 (0, 0) にするなどしてみてください．

質点間の衝突はまだ発生しない点に注意してください．

PointMass を継承した FixedPointMass クラスも用意されています．固定端につながれたバネ質量は以下のように配置できます．

simple_main.py (ver 5.0)
32
33    p3 = spm.FixedPointMass((350, 50), (0, 0), world, drawer=circle_drawer)
34    p4 = spm.PointMass((450, 80), (0, 0), world, drawer=circle_drawer)
35    sp34 = spm.Spring(p3, p4, world, spring_const, natural_len, line_drawer)
36    actor_list += [p3, p4, sp34]
37
38    while True:

7.1.4. 質点間の衝突解決

質点間の衝突を発生させるには，CollisionResolver 型のオブジェクトを actor_list に入れます．

初期化の際は，World 型オブジェクトと，衝突判定の対象を含むリストを渡します．ここに含まれる質点の相互間で衝突が発生します．

simple_main.py (ver 6.0)
37
38    actor_list.append(spm.CollisionResolver(world, actor_list))
39
40    while True:

7.1.5. 描画処理のカスタマイズ

描画処理クラスを自作する例を示します．まずは spring_mass モジュールが用意している CircleDrawer と LineDrawer の定義を見るのが早いです:

class CircleDrawer:
    def __init__(self, color, width):
        self.color = pygame.Color(color)
        self.width = width

    def __call__(self, screen, center, radius):
        pygame.draw.circle(screen, self.color, center, radius, self.width)


class LineDrawer:
    def __init__(self, color, width):
        self.color = pygame.Color(color)
        self.width = width

    def __call__(self, screen, pos1, pos2):
        pygame.draw.line(screen, self.color, pos1, pos2, self.width)

つまり，PointMass 生成時に渡す描画オブジェクトは，引数 screen, center, radius を受け取る呼び出し可能オブジェクトであり， Spring 生成時に渡す方は screen と両端の位置 pos1, pos2 を受け取る呼び出し可能オブジェクトです．

これに倣って CircleDrawer の代わりに適当な画像を表示する ImageDrawer を作ってみましょう．

simple_main.py (ver 7.0)
 1import pygame
 2import spring_mass as spm
 3
 4
 5class ImageDrawer:
 6    def __init__(self, image):
 7        self.image = image
 8
 9    def __call__(self, screen, center, radius):
10        screen.blit(self.image, (center[0] - self.image.get_width() / 2,
11                                 center[1] - self.image.get_height() / 2))
12

48
49    image = pygame.image.load("../../assets/player/p1_walk03.png").convert_alpha()
50    image_drawer = ImageDrawer(image)
51    p5 = spm.PointMass((200, 300), (5, -5), world, radius=image.get_height() / 2,
52                       mass=50, restitution=0.6, drawer=image_drawer)
53    actor_list.append(p5)
54
55    while True:

ImageDrawer の定義は 5-11行目の通りです．初期化時に画像 (Surface 型オブジェクト) を受け取り，__call__ される度に blit で貼り付けます．

49-50行目では，hello_pygame.py で使ったプレイヤ画像を読み込み，ImageDrawer 型オブジェクトを生成しています．後は CircleDrawer と同じ使い方です．

なお，画像の形がどんなものでも，衝突計算は円形として行われます．

もっといろんな形の図形の衝突計算をしたいのですが．
完全に任意の形状となると難しいです．

このテキストでは，円どうしの衝突と，円と直線の衝突のみを扱っています．これらの判定は簡単だからです．もっと複雑な形状の衝突を考えようとすると，一般論はとても難しく，それだけで本が一冊書けてしまうくらいです．

円と直線以外で比較的扱いやすいものとしては，縦横が座標軸方向と揃っている軸平行長方形 (axis-aligned rectangle, aa-rectangle) や，長方形の両端に半円をくっつけたカプセル (capsule) などがあります．ゲームなどのように厳密性が求められない場合は，複雑な形状もこれらで近似してしまうことも多いです．興味のある人は調べてみてください．

なお，pygame には軸平行長方形の衝突判定を行う sprite というモジュールが用意されています．

7.2. Slinky Drop シミュレーション

少し応用っぽいものを作ってみます．Slinky Drop という理科実験の話を聞いたことはあるでしょうか．以下の動画で紹介されています．

この動画は実験結果が判明する直前で終わります．結果が気になる人は，続きの動画を探すか，またはせっかくなので自分でシミュレーションしてみてはどうでしょう．

simple_main.py の最初のバージョン (ver 1.0) をコピーして， slinky_drop_main.py を作り，バネと質量を追加してスリンキーを模擬します．

slinky_drop_main.py (ver 8.0)
 1import pygame
 2import spring_mass as spm
 3
 4
 5def create_point_mass(pos, world, fixed=False):
 6    vel = (0, 0)
 7    radius = 10
 8    mass = 0.1
 9    viscous = 0.01
10    restitution = 0.95
11    if fixed:
12        PointMassClass = spm.FixedPointMass
13    else:
14        PointMassClass = spm.PointMass
15    return PointMassClass(pos, vel, world, radius, mass, viscous, restitution,
16                          spm.CircleDrawer("green", width=0))
17
18
19def create_spring(p1, p2, world):
20    spring_const = 0.02
21    natural_len = 5
22    return spm.Spring(p1, p2, world, spring_const, natural_len,
23                      spm.LineDrawer("white", width=2))
24
25
26def main():
27    pygame.init()
28    width, height = 600, 400
29    screen = pygame.display.set_mode((width, height))
30    clock = pygame.time.Clock()
31
32    world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.5))
33    actor_list = []
34
35    p = []
36    sp = []
37    spacing = 30
38    p.append(create_point_mass((width/2, 0), world, fixed=True))
39    for k in range(1, 15):
40        p.append(create_point_mass((width/2, spacing * k), world))
41        sp.append(create_spring(p[k - 1], p[k], world))
42    actor_list += p + sp
43
44    while True:
45        frames_per_second = 60
46        clock.tick(frames_per_second)
47
48        should_quit = False
49        for event in pygame.event.get():
50            if event.type == pygame.QUIT:
51                should_quit = True
52            elif event.type == pygame.KEYDOWN:
53                if event.key == pygame.K_ESCAPE:
54                    should_quit = True
55                elif event.key == pygame.K_SPACE:
56                    sp[0].is_alive = False
57        if should_quit:
58            break
59
60        for a in actor_list:
61            a.update()
62        actor_list[:] = [a for a in actor_list if a.is_alive]
63
64        screen.fill(pygame.Color("black"))
65        for a in actor_list:
66            a.draw(screen)
67        pygame.display.update()
68
69    pygame.quit()
70
71
72if __name__ == "__main__":
73    main()
5-16行目，19-23行目:
質点の生成，バネの生成を行う関数を用意しておきます．メインプログラム側では，これらを呼び出すことで細かい引数の設定などに煩わされることなくオブジェクトを生成できます．このようにオブジェクト生成を行う関数をファクトリ関数と呼ぶことがあります．
35-42行目
画面中央の一番上に固定質点を置き，そこから15個の質点をバネで直列につなぎます．
52-56行目
Space キーを押すと，最上部のバネの is_alive 属性が False になり，手を離したことを模擬します．

さて，予想通りの動きだったでしょうか?

56行目で sp[0].is_alive = False としているところは「属性はクラスの外から書き換えるのは避ける」というこのテキストの方針に違反してませんか?
違反ですね．

真面目に方針を守るなら，Spring クラスに自身の is_alive 属性を False にするメソッド (例えば名前は disappear とか) を設けてそれを呼ぶべきです．あるいは Actor はすべてこのメソッドを持つべきかも知れません．

「is_alive を False にするだけの処理をわざわざメソッドにする必要があるのか?」と考える人も多いと思います．実際，Python プログラマの場合はこのくらいは気にせずにクラス外から直接書き換える人も多く，その方が Python らしいという人もいます．

このテキストで「属性を外から書き換えない」方針をできるだけ守ってきたのは，オブジェクトに責務を任せる (オブジェクトの状態を変更するのはオブジェクト自身にやらせる) という考え方に慣れてもらうためです．その考え方が十分に理解できるようになったなら，実際のコードの書き方のルールは自分で(あるいはチームで) 決めればよいと思います．

7.3. インタラクティブ放物シミュレーション

もう少し複雑な例として，この章の冒頭にアニメーションを示した，マウス・キー操作でインタラクティブに質量やバネを配置できるものを作ります．

まずは空っぽのフレームワークです．ちょっと複雑になるのでメインプログラム側もクラスとして作ります．interactive_main.py として保存してください．

interactive_main.py (ver 9.0)
 1import pygame
 2import spring_mass as spm
 3
 4
 5class AppMain:
 6    def __init__(self):
 7        pygame.init()
 8        width, height = 600, 400
 9        self.screen = pygame.display.set_mode((width, height))
10
11        self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
12        self.actor_list = []
13
14    def update(self):
15        for a in self.actor_list:
16            a.update()
17        self.actor_list[:] = [a for a in self.actor_list if a.is_alive]
18
19    def draw(self):
20        self.screen.fill(pygame.Color("black"))
21        for a in self.actor_list:
22            a.draw(self.screen)
23        pygame.display.update()
24
25    def run(self):
26        clock = pygame.time.Clock()
27
28        while True:
29            frames_per_second = 60
30            clock.tick(frames_per_second)
31
32            should_quit = False
33            for event in pygame.event.get():
34                if event.type == pygame.QUIT:
35                    should_quit = True
36                elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
37                    should_quit = True
38            if should_quit:
39                break
40
41            self.update()
42            self.draw()
43
44        pygame.quit()
45
46
47if __name__ == "__main__":
48    AppMain().run()

7.3.1. 質点と壁と衝突処理

まず，質点と壁を生成し，衝突処理を発動させる部分だけ作ります．これらはオブジェクトを配置するだけなので，比較的簡単です．

interactive_main.py (ver 10.0)
 1import random
 2import pygame
 3import spring_mass as spm
 4
 5
 6class ActorFactory:
 7    def __init__(self, world, actor_list):
 8        self.world = world
 9        self.actor_list = actor_list
10
11    def create_point_mass(self, pos, fixed=False):
12        vel = (random.uniform(-10, 10), random.uniform(-10, 0))
13        mass = 10
14        radius = 10
15        viscous = 0.01
16        restitution = 0.95
17        if fixed:
18            PointMassClass = spm.FixedPointMass
19            color = "gray"
20        else:
21            PointMassClass = spm.PointMass
22            color = "green"
23        return PointMassClass(pos, vel, self.world, radius, mass, viscous,
24                              restitution, spm.CircleDrawer(color, width=0))
25
26    def create_collision_resolver(self):
27        return spm.CollisionResolver(self.world, self.actor_list)
28
29    def create_boundary(self, name):
30        width, height = self.world.size
31        geometry = {"top": ((0, -1), (0, 0)),
32                    "bottom": ((0, 1), (0, height)),
33                    "left": ((-1, 0), (0, 0)),
34                    "right": ((1, 0), (width, 0))}
35        normal, point_included = geometry[name]
36        return spm.Boundary(normal, point_included, self.world, self.actor_list)
37
38
39class AppMain:
40    def __init__(self):
41        pygame.init()
42        width, height = 600, 400
43        self.screen = pygame.display.set_mode((width, height))
44
45        self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
46        self.actor_list = []
47        self.factory = ActorFactory(self.world, self.actor_list)
48
49        self.actor_list.append(self.factory.create_collision_resolver())
50        self.actor_list.append(self.factory.create_boundary("top"))
51        self.actor_list.append(self.factory.create_boundary("bottom"))
52        self.actor_list.append(self.factory.create_boundary("left"))
53        self.actor_list.append(self.factory.create_boundary("right"))
54
55    def update(self):

73            should_quit = False
74            for event in pygame.event.get():
75                if event.type == pygame.QUIT:
76                    should_quit = True
77                elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
78                    should_quit = True
79                elif event.type == pygame.MOUSEBUTTONDOWN:
80                    self.actor_list.append(self.factory.create_point_mass(event.pos))
81            if should_quit:
82                break
6-36行目

利便性のため，ファクトリ関数群をメソッドとして持つクラス ActorFactory を作っておきます．まずこのオブジェクトに world や actor_list を保持させておき，Actor を生成するときに内部で渡してやるようにすれば，メインプログラム側でいちいちこれらを引数として渡す必要がなくなります．

29-36行目では Boundary のファクトリを定義していますが，上下左右の境界を法線ベクトルと点の座標で指定するのは煩雑なので，名前で指定できるようにしています．"top", "bottom", "left", "right" という文字列をキーとして，法線ベクトルと点の座標のタプルを返す辞書 geometry を作り，必要な値を入れておいてオブジェクト生成時に使います．

47-53行目
ファクトリオブジェクト factory を生成し，これを使って衝突処理オブジェクトと境界オブジェクトを生成し，actor_list に加えます．まあまあすっきり書けるようになったと思います．
79-80行目
マウスボタンが押されたら質点を生成し，actor_list に加えます．この時点では，どのボタンかに関わらず非固定の質点が生成されます．
18行目と21行目で，PointMassClass という変数に「クラス名」を代入していますが，こんなことしていいんですか?
はい，Python では問題ないです．

過去の Q&A で何度か言及しましたが，Python ではクラスもオブジェクトの一種として扱えて，変数に割り当てたり，関数に引数として渡したりできます．23行目のように変数を介してインスタンスを生成することもできます．

ここでは PointMass と FixedPointMass を生成するときの引数は全く同じでよいので，if と else のそれぞれで生成するより，このように書く方が簡潔です．

ファクトリ関数たちを AppMain のメソッドにしてしまえば，そもそも world や actor_list を渡す必要なくなって楽なんじゃないですか?
AppMain クラスの肥大化を防ぐために，外部に出せるものはできるだけ外部に出すことにしています．

何度も話しているように，プログラムの大規模化に立ち向かうときの鍵は，同時に考慮しなくてはならないものを減らすことです．クラスはできるだけ小さく保つ方が，同時に考慮すべきものが少なくなります．

ActorFactory が持つメソッドは，AppMain の属性を変更しないように作られています．こういうものはできるだけ外に追い出す方がよいです．

一方，AppMain の属性を更新するような処理は迂闊にクラス外に追い出せません (このテキストでは，属性はできるだけクラス外からは変更しない方針を採ってきたのを思い出してください)． ActorFactory が AppMain の actor_list 属性への参照を保持しているからといって，「そうだ，create_point_mass メソッドの中で actor_list に append してしまえば一石二鳥だ!」などとは考えない方がよいです．

よく知られている経験則として「オブジェクトを作るクラスとそのオブジェクトを使うクラスは分けておく方がよい」という設計指針が挙げられます．ファクトリという「オブジェクトを生成するだけ」の関数やクラスを作るのは，そういう指針の一環であるといえます．

複数のファクトリ関数をクラスにまとめている意義があまり感じられないのですが．
この例ではあまりピンと来ないかも知れません．

少しだけメリットが現れているのは，world や actor_list のように共通で必要となるようなものを毎回渡さずに済む点です．

もっと有難みが出るのは，生成したいオブジェクトの組み合わせをごっそり入れ替えられるようにしたい場合です．次の節で ActorFactory にバネを生成するメソッドを追加します．バネや質点のパラメータ群はこの章では 1 種類に固定していますが，目的に応じて異なるバネマスセット (例えば重くて硬いセットとか，柔らかく自然長の長いセットとか) を用意したい場合もあり得ます．パラメータを外部ファイルから読み込みたい場合などもあるでしょう．

そのような場合，それぞれのバネマスセットを ActorFactory の派生クラスとして定義しておけば，47行目で作る factory を変えるだけで簡単に切り替えられるようになります．

7.3.2. バネを追加

あとはバネを作る方法を加えれば完成です．Space キーを押している間は，生成した質点がバネで順につながれるようにします．

interactive_main.py (ver 11.0)
25
26    def create_spring(self, p1, p2):
27        spring_const = 0.01
28        natural_len = 20
29        return spm.Spring(p1, p2, self.world, spring_const, natural_len,
30                          spm.LineDrawer("white", width=2))
31
32    def create_collision_resolver(self):

45class AppMain:
46    def __init__(self):
47        pygame.init()
48        width, height = 600, 400
49        self.screen = pygame.display.set_mode((width, height))
50
51        self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
52        self.actor_list = []
53        self.factory = ActorFactory(self.world, self.actor_list)
54
55        self.actor_list.append(self.factory.create_collision_resolver())
56        self.actor_list.append(self.factory.create_boundary("top"))
57        self.actor_list.append(self.factory.create_boundary("bottom"))
58        self.actor_list.append(self.factory.create_boundary("left"))
59        self.actor_list.append(self.factory.create_boundary("right"))
60
61        self.point_mass_prev = None
62
63    def add_connected_point_mass(self, pos, button):
64        if button == 1:
65            fixed = False
66        elif button == 3:
67            fixed = True
68        else:
69            return
70        p = self.factory.create_point_mass(pos, fixed)
71        self.actor_list.append(p)
72
73        if self.point_mass_prev is not None:
74            sp = self.factory.create_spring(p, self.point_mass_prev)
75            self.actor_list.append(sp)
76        if pygame.key.get_pressed()[pygame.K_SPACE]:
77            self.point_mass_prev = p
78
79    def update(self):

 97            should_quit = False
 98            for event in pygame.event.get():
 99                if event.type == pygame.QUIT:
100                    should_quit = True
101                elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
102                    should_quit = True
103                elif event.type == pygame.KEYUP and event.key == pygame.K_SPACE:
104                    self.point_mass_prev = None
105                elif event.type == pygame.MOUSEBUTTONDOWN:
106                    self.add_connected_point_mass(event.pos, event.button)
107            if should_quit:
108                break
26-30行目
ファクトリに Spring を作るメソッドを加えておきます．
61行目
質点を順にバネでつなぐためには，一つ前に作った質点を覚えておく必要があります．ちょっと安直ですが，AppMain の属性 point_mass_prev として覚えておくことにします．初期値は None (覚えてない) です．
63-77行目

マウスボタンが押されたときの処理は，ちょっと複雑なのでメソッドに分けておくことにします．

まず64-71行目で，ボタン番号に応じて非固定または固定の質点を生成し， actor_list に入れておきます．

次の73-75行目で，一つ前の質点を覚えているなら，さっき作った質点とバネでつなぎます．

最後に，Space キーが押されているなら，さっき作った質点を point_mass_prev 属性として記憶しておきます．

103-104行目
忘れてはいけないのは，Space キーから手を離したときに point_mass_prev を None にリセットしておくことです．これをしないと何が起きるか想像できるでしょうか．この2行をコメントアウトして答え合わせしてみてください．
106行目
マウスボタンが押されたときは add_connected_point_mass に処理をさせます．
Space キーの処理が run と add_connected_point_mass に分散しているのがわかりにくいです．
はい，これはあまりよい設計ではないです．

この点に限らず，このテキスト全体を通して，キー・マウス入力処理は抽象化が足りていません．あまり凝った構造にしても学びにくいと考えたためです．

本来は，キー・マウス入力に関する処理は AppMain から独立させた上で，必要に応じてそちらから AppMain のメソッドを呼び出させたり，AppMain 側からの状態問い合わせを受け付けたりするように設計する方がよいでしょう．

if - elif - elif ... を羅列するのもあまり読みやすいコードとはいえません．イベント処理を行うクラスをイベントの種類ごとに用意して，event.type をキーとする辞書に入れておくとか，「Chain of Responsibility」と呼ばれるデザインパターンを利用するとか，改善方法はいろいろと考えられます．

7.4. spring_mass モジュールを作る

アプリケーション側の事例はここまでとして，これまで天下りで与えていた spring_mass モジュールを見ていきます．

7.4.1. 基本設計

まず基本的な構造として，Actor はメソッド update と draw を持ち， is_alive 属性を持つものとします．

update が呼ばれると，Actor は自身が力を及ぼすべき相手の receive_force メソッドを呼び出します．その際に引数として力ベクトルを渡します．力を及ぼすべき相手としては PointMass オブジェクトが想定されており，属性 pos, vel, radius, mass, viscous_damping, restitution を読み取れるとします．

PointMass は update メソッドによって，自身に及ぼす力を計算するとともに，受け取った力による自身の運動を計算します．

このルールの下で各 Actor を実装することで，メインプログラム側では:

for a in actor_list:
    a.update()
actor_list[:] = [a for a in actor_list if a.is_alive]

のように書くだけで，a の具体的な型 (PointMass か Spring かなど) に関わらずすべての Actor が協調して動きます．また:

for a in actor_list:
    a.draw(screen)

のように書くことで a の具体的な型に関わらず必要な描画が行われます．その「update」によってバネ復元力を発生するのか反発力を発生するのか運動計算をするのかを決めるは各 Actor の責務で，「draw」によって何が描かれるかも各 Actor の責務です．この性質をポリモーフィズムと呼ぶのでした．

Python だから PointMass でも Spring でも構わずに actor_list に入れられますが，Java や C++ のようにリストの型を指定しなくてはいけない言語ではそうはいかないと思います．

6章で Particle と ConfinedParticle を使い分けたときは，継承関係があったから基底クラスである Particle 型のリストを使えばよかったのですが，今回は PointMass や Spring などには継承関係がありません．こういうときはどうすればよいですか?

PointMass や Spring などの共通の親となるような基底クラスを作ればよいです．

この章で Actor と呼んでいるものは，実は概念的にはそのような共通の基底クラスを表しています．

Actor 自身のインスタンスを生成することはありません．派生クラスを作ることのみを目的として用意するものです．このようなクラスを抽象基底クラス (abstract base class, ABC) と呼びます．

質問者も言っているように Python では抽象基底クラスを定義する必要はありません．しかし，定義しておくとプログラムの意味が明確になり，型チェックなどもかけられるようになります．興味がある人は abc モジュールについて調べてみてください．

なお，ここで:

actor_list[:] = [a for a in actor_list if a.is_alive]

の [:] は必須です．ここを actor_list = ... としてしまうと， Boundary 型や CollisionResolver 型のオブジェクトに渡した actor_list とは別のものを指してしまうからです．

7.4.2. 質点

では spring_mass.py の内容をすべて削除するところから始めます．その空のファイルに以下を記載します．

spring_mass.py (ver 12.0)
 1import pygame
 2
 3
 4PgVector = pygame.math.Vector2
 5
 6
 7class World:
 8    def __init__(self, size, dt, gravity_acc):
 9        self.size = size
10        self.dt = dt
11        self.gravity_acc = PgVector(gravity_acc)
12
13
14class CircleDrawer:
15    def __init__(self, color, width):
16        self.color = pygame.Color(color)
17        self.width = width
18
19    def __call__(self, screen, center, radius):
20        pygame.draw.circle(screen, self.color, center, radius, self.width)
21
22
23class LineDrawer:
24    def __init__(self, color, width):
25        self.color = pygame.Color(color)
26        self.width = width
27
28    def __call__(self, screen, pos1, pos2):
29        pygame.draw.line(screen, self.color, pos1, pos2, self.width)
30

以降で pygame.math.Vector2 を多用するため，4行目で PgVector という別名で使えるようにしています (ちょっと不自然ですが，実はこれは次の章への布石です)．

World は説明不要だと思います．particle.py の World とは微妙に違うことにだけ注意してください．

CircleDrawer と LineDrawer は既に説明した通りです．

さらに，質点クラスを定義する以下のコードを追記します．

spring_mass.py (ver 12.0)
31
32def compute_gravity_force(mass, gravity_acc):
33    return mass * gravity_acc
34
35
36def compute_viscous_damping_force(viscous_damping, vel):
37    return -viscous_damping * vel
38
39
40def integrate_symplectic(pos, vel, force, mass, dt):
41    vel_new = vel + force / mass * dt
42    pos_new = pos + vel_new * dt
43    return pos_new, vel_new
44
45
46class PointMass:
47    def __init__(self, pos, vel, world, radius=10.0, mass=10.0,
48                 viscous_damping=0.01, restitution=0.95, drawer=None):
49        self.is_alive = True
50        self.world = world
51        self.drawer = drawer or CircleDrawer("blue", 1)
52
53        self.pos = PgVector(pos)
54        self.vel = PgVector(vel)
55        self.radius = radius
56        self.mass = mass
57        self.viscous_damping = viscous_damping
58        self.restitution = restitution
59
60        self.total_force = PgVector((0, 0))
61
62    def update(self):
63        self.generate_force()
64        self.move()
65        self.total_force = PgVector((0, 0))
66
67    def draw(self, screen):
68        self.drawer(screen, self.pos, self.radius)
69
70    def receive_force(self, force):
71        self.total_force += PgVector(force)
72
73    def generate_force(self):
74        force_g = compute_gravity_force(self.mass, self.world.gravity_acc)
75        force_v = compute_viscous_damping_force(self.viscous_damping, self.vel)
76        self.receive_force(force_g + force_v)
77
78    def move(self):
79        self.pos, self.vel = \
80            integrate_symplectic(self.pos, self.vel, self.total_force, self.mass, self.world.dt)

PointMass は，重力と粘性抵抗力を自身に及ぼす役割と，受け取った力に応じて位置と速度を更新する役割の 2 つを持っています．前者は update メソッドの中から generate_force メソッドを呼び出すことで行います (73-76行目)．

\[\begin{split}\boldsymbol{f}_\text{gravity} &= m \boldsymbol{g}\\ \boldsymbol{f}_\text{viscous} &= - d \boldsymbol{v}\end{split}\]

ただし \(m\), \(d\), \(\boldsymbol{v}\) はそれぞれ質量，粘性抵抗係数，速度です．計算した力ベクトルを，自身の receive_force メソッドを呼び出して渡します．

後者の役割のために，total_force という合力を蓄積しておく属性を用意し，最初はゼロベクトルにしておきます (60行目)． receive_force メソッドで受け取った力ベクトルはここに足しこんでいきます (71行目)． update メソッドでは，move メソッドを呼ぶことでその時点までの合力から運動を計算します (78-80行目)．加速度が力から計算されている以外は particle.py と同じです．

\[\begin{split}\boldsymbol{a} &= \frac{\boldsymbol{f}_\text{total}}{m}\\ \Delta \boldsymbol{v} &= \boldsymbol{a} \Delta t\\ \Delta \boldsymbol{x} &= \boldsymbol{v} \Delta t\end{split}\]

運動が確定したら，属性 total_force をゼロにリセットしておきます (65行目)．

このように receive_force を介して力計算と運動計算を分けるのは回りくどく見えるかもしれません．しかしこうすることにより，バネや壁など他の要素から受け取る力と統一的に扱えるようになります．

残るメソッドは draw です．初期化時に用意しておいた描画ストラテジ self.drawer を呼び出します．

self.drawer を初期化しているのは 51 行目ですが，ここの drawer or CircleDrawer("blue", 1) という書き方は見慣れない人が多いと思います．これは， drawer が None でなければその値にし， None ならば or の後ろの値にするという意味です．このように「もし値が未定ならばこの値にする」というときによく使う書き方です．

今の drawer or CircleDrawer("blue", 1) の部分がなぜそのような動作になるのかわかりません．
ちょっと難しいのですが，これは短絡 (shortcut) と呼ばれる仕様で，多くのプログラミング言語で採用されています．

まず or は真偽値に対して「または」の演算をすることを思い出してください．「x または y」は，x と y のどちらか一方が真であれば，もう一方の値をチェックすることなく真になります．

このことを利用して，Python を含む多くのプログラミング言語では， or 演算の対象を左から順に評価していき，真であるものが見つかった時点でその値を返し，それより後のものは無視します． and も同様で，左から順に見ていって偽が現れた時点でその値を返します．(C 言語の || や && も同じです)

なぜそのような仕様になっているかというと，以下のような書き方をするときに便利だからです:

if k >= 0 and k < len(array) and array[k] == some_value:
    do_something()

この例では，整数 k がシーケンス array のインデックスの範囲内であることを確認してから array[k] の値を調べています．もし and 演算子が左から見ていくと保証されていなかったとすると，インデックス範囲をチェックする前に array[k] を評価してエラーになってしまうかもしれません．

元の例では，Python の None は真偽値として評価されたときに偽になることに注意します．もし drawer が None なら後ろの CircleDrawer("blue", 1) に評価が移りこれが結果になります．もし drawer に描画オブジェクトが渡されている場合は偽になりませんので，後ろの CircleDrawer("blue", 1) は無視されて， drawer の値が結果になります．

このような書き方を自分でも使ってみようと思う人は，Python で偽と評価されるのは None だけではないことに注意してください．例えば 0，空文字列 ""，空リスト [] なども真偽値として評価すると偽になりますので， or 演算にかけてもそこで止まらず次に進んでしまいます．

ここまでの動作を確認するため， interactive_main.py から CollisionResolver 型オブジェクトと Boundary 型オブジェクトを生成している部分をいったんコメントアウトしておきます．

interactive_main.py (ver 12.0)
45class AppMain:
46    def __init__(self):
47        pygame.init()
48        width, height = 600, 400
49        self.screen = pygame.display.set_mode((width, height))
50
51        self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
52        self.actor_list = []
53        self.factory = ActorFactory(self.world, self.actor_list)
54
55        # self.actor_list.append(self.factory.create_collision_resolver())
56        # self.actor_list.append(self.factory.create_boundary("top"))
57        # self.actor_list.append(self.factory.create_boundary("bottom"))
58        # self.actor_list.append(self.factory.create_boundary("left"))
59        # self.actor_list.append(self.factory.create_boundary("right"))
60
61        self.point_mass_prev = None
62

これで，左クリックで質点を生成することだけはできるようになったはずです．重力加速度や初速度はお好みで調整してください．

右クリック (固定質点の生成) やスペースを押しながらのクリック (バネ生成) をするとエラーになるので注意してください．まだ作ってないからです．

7.4.3. 固定質点

次に，PointMass を継承して FixedPointMass を作ります．以下のコードを spring_mass.py に追加します．

spring_mass.py (ver 13.0)
82
83class FixedPointMass(PointMass):
84    def __init__(self, pos, vel, world, radius=10, mass=10,
85                 viscous_damping=0.01, restitution=0.95, drawer=None):
86        super().__init__(pos, vel, world, radius, mass,
87                         viscous_damping, restitution, drawer)
88        self.vel, self.mass = PgVector((0, 0)), 1e9
89
90    def move(self):
91        pass

move メソッドを上書きして動かなくする (90-91行目) だけで済みそうに思うのですが，実は速度と質量も上書きしておかないと衝突計算が正しくできません．

PointMass の __init__ メソッドの内容をコピーペーストするのは無駄が多いので，PointMass の __init__ メソッドを呼び出した上で， self.vel にゼロベクトルを， self.mass に非常に大きな値を代入し直すことにします．

super().__init__ は継承元 (スーパークラス) の __init__ メソッドを表す Python の記法です．

7.4.4. バネ

続いてバネを定義します．

spring_mass.py (ver 14.0)
 93
 94def compute_restoring_force(pos1, pos2, spring_const, natural_len):
 95    if pos1 == pos2:
 96        return None
 97    vector12 = pos2 - pos1
 98    distance = vector12.magnitude()
 99    unit_vector12 = vector12 / distance
100    f1 = unit_vector12 * spring_const * (distance - natural_len)
101    return f1
102
103
104class Spring:
105    def __init__(self, point_mass1, point_mass2, world,
106                 spring_const=0.01, natural_len=0.0, drawer=None):
107        self.is_alive = True
108        self.world = world
109        self.drawer = drawer or LineDrawer("blue", 1)
110
111        self.p1 = point_mass1
112        self.p2 = point_mass2
113        self.spring_const = spring_const
114        self.natural_len = natural_len
115
116    def update(self):
117        if not (self.p1.is_alive and self.p2.is_alive):
118            self.is_alive = False
119            return
120        self.generate_force()
121
122    def draw(self, screen):
123        self.drawer(screen, self.p1.pos, self.p2.pos)
124
125    def generate_force(self):
126        f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
127        if f1 is None:
128            return
129        self.p1.receive_force(f1)
130        self.p2.receive_force(-f1)

Spring クラスの update メソッドでは，まず両端の質点の存在を確認して，どちらかが存在しなければ自身の is_alive を False にします (117-119行目)．両方存在するなら generate_force メソッドで力を及ぼします．

バネ復元力の計算自体は難しいことはないです．質点1が受け取る力は

\[\boldsymbol{f}_1 = k (L - L_\text{natural}) \boldsymbol{n}_\text{1,2}\]

です．ただし \(\boldsymbol{n}_\text{1,2}\) は質点1から2に向かう単位ベクトルで，\(k\), \(L\), \(L_\text{natural}\) はそれぞれバネ定数，現在の長さ，自然長です．

95-100行目は，基本的にはこの式の通りですが，注意が必要なのは単位ベクトルを求めるために距離で割っているところで (99行目)，距離が 0 だとゼロ除算が発生します．質点が有限の大きさを持っていて (だからそれは質点じゃないだろというツッコミはスルーします) 衝突判定が有効なら本来この状況は生じませんが，シミュレーション上では生じ得ます．ここでは単純に，その場合は力を発生しないことにします．

interactive_main.py を実行して，Space キーを押しながら質点を配置してバネの動作を確認してください．重力加速度や初速度をゼロにしておく方がわかりやすいかも知れません．

単体テストはどういうふうに書けばよいですか．
興味がある人はつづきをどうぞ (test_spring_mass.py ver 14.0)．

test_spring_mass.py を作成して以下を記載し，Explorer サイドバーから Run All Tests を実行してください．このフォルダではまだ pytest が有効化されていないので，有効化を促すポップアップが出るはずです．particle フォルダと同じく有効化してください．

関数 generate_restoring_force をテストする関数を作ります．バネ復元力は「正解の値」を用意するのが簡単なので，引数 f1_expected として与えてやることにしています．また，None が返ることが期待される場合 (復元力の計算しない場合) は引数 expects_none を True にします．これらの引数は @pytest.mark.parametrize で与えます．

x 方向と y 方向それぞれに伸びているバネしかテストしていません．本当は斜め方向もテストするべきですが，被テスト関数側ではベクトル表記で計算しているので，独立な 2 方向が正しく動けばまあ大丈夫かなということで省略しています．気になる人は，45度や60度などのわかりやすい方向のテストを加えてみるとよいと思います．

test_spring_mass.py (ver 14.0)
 1import pytest
 2import spring_mass as spm
 3from spring_mass import PgVector
 4
 5
 6@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
 7    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
 8    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
 9    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
10    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),
11    ((0, 0), (10, 0), 2, 13,  (-2 * (13 - 10), 0), False),
12    ((10, 0), (0, 0), 2, 13,  (2 * (13 - 10), 0), False),
13
14    ((0, 0), (0, 10), 3, 0,  (0, 3 * 10), False),
15    ((2, 10), (2, 5), 4, 0,  (0, -5 * 4), False),
16    ((0, 0), (0, 10), 3, 5,  (0, 3 * 5), False),
17    ((0, 0), (0, 10), 3, 12,  (0, -3 * 2), False),
18
19    ((0, 0), (0, 0), 0, 0, (0, 0), True), # equal positions
20    ((12, 5), (12, 5), 0, 0, (0, 0), True), # equal positions
21])
22def test_restoring_force(pos1v, pos2v, k, nl, f1_expected, expects_none):
23    pos1 = PgVector(pos1v)
24    pos2 = PgVector(pos2v)
25    f1 = spm.compute_restoring_force(pos1, pos2, k, nl)
26    if expects_none:
27        assert f1 is None
28    else:
29        assert f1 == PgVector(f1_expected)
こうやって数式が与えられれば，ああ確かにこういうプログラムになるなと理解できるのですが，自分で数式から考えろと言われると難しいです．特に，自然長を境に符号をどうすればよいかなど，自分で考えられる自信がありません．
単体テストを利用して，テストしながら少しずつ作っていくとよいかも知れません．

意外に聞こえるかもしれませんが，まずテストを書き，テストに通るようにプログラム本体を作っていくというやり方が有効にはたらく場合があります．テスト駆動開発 (test driven development, TDD) と呼ばれます．

バネ復元力の例で考えてみましょう．一つ前の Q&A の単体テストの説明を読んでおいてください．

さて，自然長が 0 の場合であれば，特に難しいことは考えずに作れると思います．だからまずその場合に限定したテストのみを用意します． x 方向に伸びる場合だけを考えましょう．期待される答え f1_expected は暗算で出せますね:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
])
def test_restoring_force(pos1v, pos2v, k, nl, f1_expected, expects_none):
    pos1 = PgVector(pos1v)
    pos2 = PgVector(pos2v)
    f1 = spm.compute_restoring_force(pos1, pos2, k, nl)
    if expects_none:
        assert f1 is None
    else:
        assert f1 == PgVector(f1_expected)

そして，このテストを通すような compute_restoring_force を作るのも簡単です:

def compute_restoring_force(pos1, pos2, spring_const, natural_len):
    vector12 = pos2 - pos1
    distance = vector12.magnitude()
    unit_vector12 = vector12 / distance
    f1 = unit_vector12 * spring_const * distance
    return f1

次に自然長が 0 でない場合のテストを追加していきます:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),

バネ長 = 10 よりも自然長 nl が短い場合を用意しました．これらの場合に発生する力の符号がどちらになるかは，簡単にわかるので f1_expected として指定しておきます．

先ほど作った compute_restoring_force では新しい 2 ケースのテストは通らないはずです．通るように直すのは難しくありません:

def compute_restoring_force(pos1, pos2, spring_const, natural_len):
    vector12 = pos2 - pos1
    distance = vector12.magnitude()
    unit_vector12 = vector12 / distance
    f1 = unit_vector12 * spring_const * (distance - natural_len)
    return f1

同様に，バネ長より自然長が長い場合のテストを追加します:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),
    ((0, 0), (10, 0), 2, 13,  (-2 * (13 - 10), 0), False),
    ((10, 0), (0, 0), 2, 13,  (2 * (13 - 10), 0), False),

テストを走らせてみると，これら新しいテストケースも通るはずです．さっきのコードのままでよさそうです．

y 方向のテストも加えてみます．ついでにちょっとずれた位置なんかも試してみましょう:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),
    ((0, 0), (10, 0), 2, 13,  (-2 * (13 - 10), 0), False),
    ((10, 0), (0, 0), 2, 13,  (2 * (13 - 10), 0), False),

    ((0, 0), (0, 10), 3, 0,  (0, 3 * 10), False),
    ((2, 10), (2, 5), 4, 0,  (0, -5 * 4), False),
    ((0, 0), (0, 10), 3, 5,  (0, 3 * 5), False),
    ((0, 0), (0, 10), 3, 12,  (0, -3 * 2), False),

特に問題なく通るようです．

ここまで来ると，基本的なロジックは大丈夫そうだなと自信が持てるようになります．気になる人は斜め方向のテストを加えてもよいでしょう．

あとは，例外的な条件が存在しないかを確認します．作成した compute_restoring_force を読み直すと，distance で割っているところがあるので，ゼロ除算のおそれがあります．テストケースとして追加しましょう:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),
    ((0, 0), (10, 0), 2, 13,  (-2 * (13 - 10), 0), False),
    ((10, 0), (0, 0), 2, 13,  (2 * (13 - 10), 0), False),

    ((0, 0), (0, 10), 3, 0,  (0, 3 * 10), False),
    ((2, 10), (2, 5), 4, 0,  (0, -5 * 4), False),
    ((0, 0), (0, 10), 3, 5,  (0, 3 * 5), False),
    ((0, 0), (0, 10), 3, 12,  (0, -3 * 2), False),

    ((0, 0), (0, 0), 0, 0, (0, 0), True), # equal positions
    ((12, 5), (12, 5), 0, 0, (0, 0), True), # equal positions

これらのテストケースは通りません．距離が 0 のときは None を返すように修正します:

def compute_restoring_force(pos1, pos2, spring_const, natural_len):
    if pos1 == pos2:
        return None
    vector12 = pos2 - pos1
    distance = vector12.magnitude()
    unit_vector12 = vector12 / distance
    f1 = unit_vector12 * spring_const * (distance - natural_len)
    return f1

これで ver 14.0 のコードにたどり着きました．

コツは，まず最初は例外的なケースや細かい検討を要するケースを無視して，問題のコアな部分を素直に実装することです．その後，細かいケースのテストを追加しながら精密化していくようにします．

7.4.5. 破断するバネ

Spring を拡張する例として，あるしきい値を超える力がかかると破断するバネ FragileSpring を作ってみます．

spring_mass.py (ver 15.0)
132
133class FragileSpring(Spring):
134    def __init__(self, point_mass1, point_mass2, world,
135                 spring_const=0.01, natural_len=0.0, drawer=None,
136                 break_threshold=1e9):
137        super().__init__(point_mass1, point_mass2, world, spring_const,
138                         natural_len, drawer)
139        self.break_threshold = break_threshold
140
141    def generate_force(self):
142        f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
143        if f1 is None:
144            return
145        self.p1.receive_force(f1)
146        self.p2.receive_force(-f1)
147        if f1.magnitude() > self.break_threshold:
148            self.is_alive = False

FragileSpring の __init__ メソッドでは，スーパークラスの __init__ を呼んだ上で，break_threshold 属性を追加しています (137-139行目)．

generate_force メソッドでは，スーパークラスの generate_force を呼び出して…としたいところなのですが，力の計算結果を得る手段がないので，仕方なく同じ処理内容をコピーペーストして作ります．その上で破断処理 (147-148行目) を追加します．

interactive_main.py で，ファクトリの create_spring メソッドが FragileSpring を生成するようにしておきましょう．

interactive_main.py (ver 15.0)
26    def create_spring(self, p1, p2):
27        spring_const = 0.01
28        natural_len = 20
29        break_threshold = 5.0
30        return spm.FragileSpring(p1, p2, self.world, spring_const, natural_len,
31                                 spm.LineDrawer("white", width=2), break_threshold)

このくらいの定数設定だと，画面の端から端までの長さのバネを作ると消えてなくなるはずです．

7.4.6. 衝突解決

次は CollisionResolver です．難しそうに感じますが，基本は他の Actor と同じで，力を及ぼすべき相手の receive_force メソッドを呼んで力ベクトルを渡していくのが責務です．

spring_mass.py (ver 16.0)
150
151def is_point_mass(actor):
152    return isinstance(actor, PointMass)
153
154
155def compute_impact_force_between_points(p1, p2, dt):
156    if (p1.pos - p2.pos).magnitude() > p1.radius + p2.radius:
157        return None
158    if p1.pos == p2.pos:
159        return None
160    normal = (p2.pos - p1.pos).normalize()
161    v1 = p1.vel.dot(normal)
162    v2 = p2.vel.dot(normal)
163    if v1 < v2:
164        return None
165    e = p1.restitution * p2.restitution
166    m1, m2 = p1.mass, p2.mass
167    f1 = normal * (-(e + 1) * v1 + (e + 1) * v2) / (1/m1 + 1/m2) / dt
168    return f1
169
170
171class CollisionResolver:
172    def __init__(self, world, actor_list, target_condition=None, drawer=None):
173        self.is_alive = True
174        self.world = world
175        self.drawer = drawer or (lambda surface: None)
176
177        self.actor_list = actor_list
178        if target_condition is None:
179            self.target_condition = is_point_mass
180        else:
181            self.target_condition = target_condition
182
183    def update(self):
184        self.generate_force()
185
186    def draw(self, surface):
187        self.drawer(surface)
188
189
190    def generate_force(self):
191        plist = [a for a in self.actor_list if self.target_condition(a)]
192        n = len(plist)
193        for i in range(n):
194            for j in range(i + 1, n):
195                p1, p2 = plist[i], plist[j]
196                f1 = compute_impact_force_between_points(p1, p2, self.world.dt)
197                if f1 is None:
198                    continue
199                p1.receive_force(f1)
200                p2.receive_force(-f1)

175 行目で描画ストラテジ self.drawer を初期化しています．引数 drawer が None でなければその値を， None なら「何も描画しない関数」をラムダ式で指定しています．

update メソッドは generate_force メソッドを呼び出すだけです (184行目)．generate_force メソッドは，まず力を及ぼすべき相手をリストアップし (191行目)，それらの相互間の衝突判定をし，衝突しているなら力を発生します．

191行目では，actor_list 属性に入っている Actor のうち対象となる条件を満たすものだけを変数 plist に取り出します．条件は target_condition 関数で判定しますが，デフォルト (178-179行目) では PointMass クラスのインスタンスかどうか (151-152行目) で判断します．条件は Strategy パターンで与えているので必要に応じて入れ替えられます．

ここで使っている isinstance(x, A) は Python が用意する関数で，オブジェクト x がクラス A のインスタンスなら True になります．ただし，A 自体ではなく A のサブクラスのインスタンスでも True になるのが特徴です．なので，FixedPointMass も衝突計算の対象になります．

plist に対象オブジェクト \(n\) 個がリストアップできたら，そのうち異なる 2 個を取り出す組合せを列挙します．一方を plist の \(i\) 番目，他方を \(j\) 番目とすると， \(i\) として \(0\) から \(n-1\) までを列挙しつつ， \(j\) として \(i + 1\) から \(n-1\) まで列挙することで重複の無い全組合せが得られます (193-194行目)．

衝突計算は compute_impact_force_between_points 関数で行います (155-168行目)．まず，中心間の距離が半径の和より大きければ衝突はしていませんので除外します (156-157行目)．衝突している場合は，中心間を結ぶベクトルに沿った方向に弾性衝突による撃力が発生するとします．

さて，ここからしばらく高校物理の時間です．質点1から質点2に向かう方向の成分だけの 1 次元運動を考えます．質点1と質点2の質量と速度をそれぞれ \(m_1\), \(m_2\), \(v_1\), \(v_2\) とすると，運動量保存則は

\[m_1 v_1 + m_2 v_2 = m_1 (v_1 + \Delta v_1) + m_2 (v_2 + \Delta v_2)\]

と表され，弾性反発係数を \(e\) とすると

\[(v_1 + \Delta v_1) - (v_2 + \Delta v_2) = - e (v_1 - v_2)\]

を満たす必要があります．そしてこれらを満たす速度変化を引き起こすために質点1と2に加えられるべき力は以下で与えられます．

\[\begin{split}f_1 &= m_1 \frac{\Delta v_1}{\Delta t}\\ f_2 &= - f_1\end{split}\]

これらを解くと以下が得られます．

\[f_1 = \frac{-(e + 1) v_1 + (e + 1) v_2}{(\frac{1}{m_1} + \frac{1}{m_2}) \Delta t}\]

というわけで，この式に従って，質点1から2に向かう方向の力を発生することになります．この方向の単位ベクトルを求め (160 行目)，その方向の速度成分を pygame.math.Vector2 の dot メソッドによる内積計算で求めます (161-162行目)．あとはこの単位方向ベクトルに力の大きさを乗じれば OK です (167行目)．

ただし，いくつか考慮すべき細かい点があります．

  • まず，質点1と2の位置が完全に一致していたら力の方向が定まりません．バネのときと同様にこれも本来はあり得ない状況ですが，計算上は起こり得ます．やはり単に無視することにします (158-159行目)．
  • 次に，particle.py で壁との衝突を考えたときにも発生した，衝突判定が複数回生じてしまう現象への対策も必要です．particle.py と同じく速度条件を考えます．つまり両者が離れていく相対速度になっているときは除外します (163-164行目)．
  • 最後に，反発係数をどのように与えるかです．質点1と質点2のそれぞれに属性 restitution が設けられていますが，本来，反発係数は衝突する二体の組ごとに与えられるもので，質点単体に定義できるものではありません．しかし，衝突の可能性のある二体ごとに利用者が定義するのは面倒なので，多くの物理シミュレータでは，質点単体に反発係数 (っぽい何か) を定義し，そこから何らかのルールで計算に用いる反発係数を求めるという便宜的方法を採用しています．ここでは両者の restitution 属性の積を使うことにします．物理的根拠はありません．

ここまでできたら，interactive_main.py に CollisionResolver を復帰させて動作を確認してください．

interactive_main.py (ver 16.0)
54        self.factory = ActorFactory(self.world, self.actor_list)
55
56        self.actor_list.append(self.factory.create_collision_resolver())
57        # self.actor_list.append(self.factory.create_boundary("top"))
58        # self.actor_list.append(self.factory.create_boundary("bottom"))
バネのときも思ったのですが，力を発生しないときに None を返すのはなぜですか? PgVector((0, 0)) を返すようにすればもっとすっきりしたコードになりますよね?
遅いからです．

バネの場合はあまり問題にならないのですが，衝突の場合はチェックする組み合わせの数が多くなります (質点の数の2乗に比例します)．そして通常ほとんどの場合に衝突は起きておらず，本来何もする必要がありません．PgVector((0, 0)) を返しても計算結果は同じですが， receive_force メソッドの呼び出しが無駄に頻繁に発生して遅くなります．

テストはどう書いたらいいですかね?
つづきをどうぞ (test_spring_mass.py ver 16.0)．

バネ復元力より数式が複雑なので，個々の「正解」を外から与えてやるのは面倒です．かといって，被テスト関数 compute_impact_force_between_points の中で計算している式をそのままテスト関数に書いても，テストの意義が弱いです．

ここではちょっと発想を変えて，被テスト関数が返す力によって変化する前後の速度を計算し，運動量が保存されていることと，弾性衝突の関係が成り立っていることを確認することにします．

test_spring_mass.py (ver 16.0)
31
32@pytest.mark.parametrize("p1v, p2v, expects_none", [
33    (((0, 0), (5, 0), 20, 10, 1.0),   # (pos, vel, r, m, e)
34     ((20, 0), (-5, 0), 20, 10, 1.0), False),
35    (((0, 0), (5, 0), 20, 15, 0.2),
36     ((20, 0), (-5, 0), 20, 30, 0.8), False),
37
38    (((0, 0), (5, 0), 20, 10, 1.0),
39     ((120, 0), (-5, 0), 20, 10, 1.0), True), # not in contact
40    (((0, 0), (5, 0), 20, 10, 1.0),
41     ((0, 0), (-5, 0), 20, 10, 1.0), True), # equal positions
42    (((0, 0), (-5, 0), 20, 10, 1.0),
43     ((20, 0), (5, 0), 20, 10, 1.0), True), # moving apart
44])
45def test_impact_force_between_points(p1v, p2v, expects_none):
46    dt, g = 1.0, 0
47    w = spm.World((600, 400), dt, g)
48    pos1, vel1, r1, m1, e1 = p1v
49    pos2, vel2, r2, m2, e2 = p2v
50
51    p1 = spm.PointMass(pos1, vel1, w,
52                radius=r1, mass=m1, restitution=e1)
53    p2 = spm.PointMass(pos2, vel2, w,
54                radius=r2, mass=m2, restitution=e2)
55    f1 = spm.compute_impact_force_between_points(p1, p2, dt)
56    if expects_none:
57        assert f1 is None
58        return
59
60    vel1_old = PgVector(vel1)
61    vel2_old = PgVector(vel2)
62    vel1_new = vel1_old + f1 / m1 * dt
63    vel2_new = vel2_old - f1 / m2 * dt
64    e = e1 * e2
65    assert vel1_new - vel2_new == -e * (vel1_old - vel2_old)
66    assert m1 * vel1_new + m2 * vel2_new == m1 * vel1_old + m2 * vel2_old

7.4.7. 壁・床

そろそろ疲れてきたと思いますが，残りは Boundary クラスです．

spring_mass.py (ver 17.0)
202
203def compute_impact_force_by_fixture(p, normal, point_included, dt):
204    invasion = normal.dot(p.pos - point_included)
205    if invasion + p.radius > 0 and normal.dot(p.vel) > 0:
206        e = p.restitution
207        v = normal.dot(p.vel)
208        m = p.mass
209        f = normal * (-(e + 1) * v) * m / dt
210    else:
211        f = None
212    return f
213
214
215class Boundary:
216    def __init__(self, normal, point_included, world, actor_list,
217                 target_condition=None, drawer=None):
218        self.is_alive = True
219        self.world = world
220        self.drawer = drawer or (lambda surface: None)
221
222        self.normal = PgVector(normal).normalize()
223        self.point_included = PgVector(point_included)
224        self.actor_list = actor_list
225        if target_condition is None:
226            self.target_condition = is_point_mass
227        else:
228            self.target_condition = target_condition
229
230    def update(self):
231        self.generate_force()
232
233    def draw(self, surface):
234        self.drawer(surface)
235
236
237    def generate_force(self):
238        plist = [a for a in self.actor_list if self.target_condition(a)]
239        for p in plist:
240            f = compute_impact_force_by_fixture(p, self.normal, self.point_included, self.world.dt)
241            if f is None:
242                continue
243            p.receive_force(f)

先ほどの CollisionResolver と似た構造をしていますが少し簡単です． generate_force メソッドでは衝突対象を列挙し (238行目)，それら1つずつに対して衝突計算をします．

境界面の正規化した法線ベクトルを \(\boldsymbol{n}\)，境界面に含まれる 1 点の座標を \(\boldsymbol{p}_0\) とします．法線ベクトルは壁面の内部に向かう方向に取っていたことを思い出しておいてください．

このとき，質点の中心座標 \(\boldsymbol{x}\) が壁内への侵入量の法線方向成分は内積 \(\boldsymbol{n} \cdot (\boldsymbol{x} - \boldsymbol{p}_0)\) により求められます．これに半径分のオフセットを加えたものが正であり，かつ速度の法線方向成分が正のとき，衝突が生じたと判定します (204-205行目)．

その際の撃力は，改めて計算してもよいのですが，質点間の衝突のときに求めた式で，相手方の速度を \(v_2 = 0\)，質量を無限大とした極限 \(m_2 \rightarrow \infty\) を考えるのが手っ取り早いです (206-209行目)．

interactive_main.py で，壁や床を作ったり外したりして遊びましょう．

interactive_main.py (ver 17.0)
56        self.actor_list.append(self.factory.create_collision_resolver())
57        self.actor_list.append(self.factory.create_boundary("top"))
58        self.actor_list.append(self.factory.create_boundary("bottom"))
59        self.actor_list.append(self.factory.create_boundary("left"))
60        self.actor_list.append(self.factory.create_boundary("right"))
61
テストは?
つづきをどうぞ (test_spring_mass.py ver 17.0)．

こちらも被テスト関数が返す力をそのままテストしても面白くないので，力によって変化する前後の速度を計算し，弾性衝突の関係が成り立っていることを確認することにします．

test_spring_mass.py (ver 17.0)
68
69@pytest.mark.parametrize("pv, n, pinc, vel_expected, expects_none", [
70    (((600, 200), (5, 3), 20, 10, 1.0),   # (pos, vel, r, m, e)
71     (1, 0), (600, 0), (-5, 3), False),
72    (((15, 200), (-5, 3), 20, 10, 0.8),
73     (-1, 0), (0, 0), (0.8 * 5, 3), False),
74    (((300, 400), (-5, 3), 20, 10, 0.6),
75     (0, 1), (0, 400), (-5, - 0.6 * 3), False),
76
77    (((600, 200), (-5, 3), 20, 10, 1.0),
78     (1, 0), (600, 0), (0, 0), True), # not moving toward warll
79    (((580, 200), (5, 3), 20, 10, 1.0),
80     (1, 0), (600, 0), (0, 0), True), # not in contact
81])
82def test_impact_force_by_fixture(pv, n, pinc, vel_expected, expects_none):
83    dt, g = 1.0, 0
84    w = spm.World((600, 400), dt, g)
85    pos, vel, r, m, e = pv
86    p = spm.PointMass(pos, vel, w, radius=r, mass=m, restitution=e)
87    normal = PgVector(n)
88    point_included = PgVector(pinc)
89    f = spm.compute_impact_force_by_fixture(p, normal, point_included, dt)
90    if expects_none:
91        assert f is None
92        return
93
94    vel_old = PgVector(vel)
95    vel_new = vel_old + f / m * dt
96    assert vel_new == PgVector(vel_expected)

7.4.8. 床面への潜り込みの回避

最後に一つ問題が残っています．particle.py でも問題になった，床面への沈み込みです．

particle.py のときと同様に y 座標を強制的に床面にすることで回避したいと思いますが，今回は少し微妙な問題が生じます．これは PointMass の責務でしょうか．Boundary の責務でしょうか．

質点の位置計算に関することなので，本来的には PointMass の責務です．しかし，PointMass 型オブジェクトは自身が床と衝突したことを知る方法がありません．バネからの力も，他の質点からの反発力も，床からの反発力も，区別せずに receive_force メソッドで受け取っているだけだからです．

では Boundary が質点の y 座標を変更するべきかというと，オブジェクトの属性を外から変更するのは避けるという方針に反するのであまりやりたくありません．

spring_mass モジュールをさらに拡充していこうと思うと，同じような問題は他の場合にも起こり得ます．例えば，他の質点と一定回数衝突したら壊れる FragilePointMass を作りたいとしましょう．しかし他の質点との衝突によって力を受け取ったことを，それ以外の力の発生と区別する手段がありません．

要するに，generate_force する側から receive_force する側に，力以外の情報を伝達する手段が必要そうです．ここでは以下のような仕組みを作ります．

spring_mass.py (ver 18.0)
  1import pygame
  2
  3
  4PgVector = pygame.math.Vector2
  5
  6
  7class World:
  8    def __init__(self, size, dt, gravity_acc):
  9        self.size = size
 10        self.dt = dt
 11        self.gravity_acc = PgVector(gravity_acc)
 12
 13
 14class CircleDrawer:
 15    def __init__(self, color, width):
 16        self.color = pygame.Color(color)
 17        self.width = width
 18
 19    def __call__(self, screen, center, radius):
 20        pygame.draw.circle(screen, self.color, center, radius, self.width)
 21
 22
 23class LineDrawer:
 24    def __init__(self, color, width):
 25        self.color = pygame.Color(color)
 26        self.width = width
 27
 28    def __call__(self, screen, pos1, pos2):
 29        pygame.draw.line(screen, self.color, pos1, pos2, self.width)
 30
 31
 32def compute_gravity_force(mass, gravity_acc):
 33    return mass * gravity_acc
 34
 35
 36def compute_viscous_damping_force(viscous_damping, vel):
 37    return -viscous_damping * vel
 38
 39
 40def integrate_symplectic(pos, vel, force, mass, dt):
 41    vel_new = vel + force / mass * dt
 42    pos_new = pos + vel_new * dt
 43    return pos_new, vel_new
 44
 45
 46class PointMass:
 47    def __init__(self, pos, vel, world, radius=10.0, mass=10.0,
 48                 viscous_damping=0.01, restitution=0.95, drawer=None):
 49        self.is_alive = True
 50        self.world = world
 51        self.drawer = drawer or CircleDrawer("blue", 1)
 52
 53        self.pos = PgVector(pos)
 54        self.vel = PgVector(vel)
 55        self.radius = radius
 56        self.mass = mass
 57        self.viscous_damping = viscous_damping
 58        self.restitution = restitution
 59
 60        self.total_force = PgVector((0, 0))
 61        self.message_list = []
 62
 63    def update(self):
 64        self.generate_force()
 65        self.move()
 66        self.total_force = PgVector((0, 0))
 67        self.message_list.clear()
 68
 69    def draw(self, screen):
 70        self.drawer(screen, self.pos, self.radius)
 71
 72    def receive_force(self, force):
 73        self.total_force += PgVector(force)
 74
 75    def receive_message(self, msg):
 76        self.message_list.append(msg)
 77
 78    def generate_force(self):
 79        force_g = compute_gravity_force(self.mass, self.world.gravity_acc)
 80        force_v = compute_viscous_damping_force(self.viscous_damping, self.vel)
 81        self.receive_force(force_g + force_v)
 82
 83    def move(self):
 84        self.pos, self.vel = \
 85            integrate_symplectic(self.pos, self.vel, self.total_force, self.mass, self.world.dt)
 86
 87        for msg in self.message_list:
 88            if msg["type"] == "floor_hit" and self.vel.y > 0:
 89                # constrain y on or above floor
 90                self.pos.y = msg["y"] - self.radius
 91
 92
 93class FixedPointMass(PointMass):
 94    def __init__(self, pos, vel, world, radius=10.0, mass=10.0,
 95                 viscous_damping=0.01, restitution=0.95, drawer=None):
 96        super().__init__(pos, vel, world, radius, mass,
 97                         viscous_damping, restitution, drawer)
 98        self.vel, self.mass = PgVector((0, 0)), 1e9
 99
100    def move(self):
101        pass
102
103
104def compute_restoring_force(pos1, pos2, spring_const, natural_len):
105    if pos1 == pos2:
106        return None
107    vector12 = pos2 - pos1
108    distance = vector12.magnitude()
109    unit_vector12 = vector12 / distance
110    f1 = unit_vector12 * spring_const * (distance - natural_len)
111    return f1
112
113
114class Spring:
115    def __init__(self, point_mass1, point_mass2, world,
116                 spring_const=0.01, natural_len=0.0, drawer=None):
117        self.is_alive = True
118        self.world = world
119        self.drawer = drawer or LineDrawer("blue", 1)
120
121        self.p1 = point_mass1
122        self.p2 = point_mass2
123        self.spring_const = spring_const
124        self.natural_len = natural_len
125
126    def update(self):
127        if not (self.p1.is_alive and self.p2.is_alive):
128            self.is_alive = False
129            return
130        self.generate_force()
131
132    def draw(self, screen):
133        self.drawer(screen, self.p1.pos, self.p2.pos)
134
135    def generate_force(self):
136        f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
137        if f1 is None:
138            return
139        self.p1.receive_force(f1)
140        self.p2.receive_force(-f1)
141
142
143class FragileSpring(Spring):
144    def __init__(self, point_mass1, point_mass2, world,
145                 spring_const=0.01, natural_len=0.0, drawer=None,
146                 break_threshold=1e9):
147        super().__init__(point_mass1, point_mass2, world, spring_const,
148                         natural_len, drawer)
149        self.break_threshold = break_threshold
150
151    def generate_force(self):
152        f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
153        if f1 is None:
154            return
155        self.p1.receive_force(f1)
156        self.p2.receive_force(-f1)
157        if f1.magnitude() > self.break_threshold:
158            self.is_alive = False
159
160
161def is_point_mass(actor):
162    return isinstance(actor, PointMass)
163
164
165def compute_impact_force_between_points(p1, p2, dt):
166    if (p1.pos - p2.pos).magnitude() > p1.radius + p2.radius:
167        return None
168    if p1.pos == p2.pos:
169        return None
170    normal = (p2.pos - p1.pos).normalize()
171    v1 = p1.vel.dot(normal)
172    v2 = p2.vel.dot(normal)
173    if v1 < v2:
174        return None
175    e = p1.restitution * p2.restitution
176    m1, m2 = p1.mass, p2.mass
177    f1 = normal * (-(e + 1) * v1 + (e + 1) * v2) / (1/m1 + 1/m2) / dt
178    return f1
179
180
181class CollisionResolver:
182    def __init__(self, world, actor_list, target_condition=None, drawer=None):
183        self.is_alive = True
184        self.world = world
185        self.drawer = drawer or (lambda surface: None)
186
187        self.actor_list = actor_list
188        if target_condition is None:
189            self.target_condition = is_point_mass
190        else:
191            self.target_condition = target_condition
192
193    def update(self):
194        self.generate_force()
195
196    def draw(self, surface):
197        self.drawer(surface)
198
199
200    def generate_force(self):
201        plist = [a for a in self.actor_list if self.target_condition(a)]
202        n = len(plist)
203        for i in range(n):
204            for j in range(i + 1, n):
205                p1, p2 = plist[i], plist[j]
206                f1 = compute_impact_force_between_points(p1, p2, self.world.dt)
207                if f1 is None:
208                    continue
209                p1.receive_force(f1)
210                p2.receive_force(-f1)
211
212
213def compute_impact_force_by_fixture(p, normal, point_included, dt):
214    invasion = normal.dot(p.pos - point_included)
215    if invasion + p.radius > 0 and normal.dot(p.vel) > 0:
216        e = p.restitution
217        v = normal.dot(p.vel)
218        m = p.mass
219        f = normal * (-(e + 1) * v) * m / dt
220    else:
221        f = None
222    return f
223
224
225class Boundary:
226    def __init__(self, normal, point_included, world, actor_list,
227                 target_condition=None, drawer=None):
228        self.is_alive = True
229        self.world = world
230        self.drawer = drawer or (lambda surface: None)
231
232        self.normal = PgVector(normal).normalize()
233        self.point_included = PgVector(point_included)
234        self.actor_list = actor_list
235        if target_condition is None:
236            self.target_condition = is_point_mass
237        else:
238            self.target_condition = target_condition
239
240    def update(self):
241        self.generate_force()
242
243    def draw(self, surface):
244        self.drawer(surface)
245
246
247    def is_floor(self):
248        return self.normal == PgVector((0, 1))
249
250    def generate_force(self):
251        plist = [a for a in self.actor_list if self.target_condition(a)]
252        for p in plist:
253            f = compute_impact_force_by_fixture(p, self.normal, self.point_included, self.world.dt)
254            if f is None:
255                continue
256            p.receive_force(f)
257            if self.is_floor():
258                p.receive_message({"type": "floor_hit", "y": self.point_included.y})
61行目, 67行目
力を受け取る側の属性として，合力とは別にメッセージを保存するリスト message_list を用意します．このリストは最初は空に初期化され (61行目)，update されるたびに空にリセットされます (67行目)．
75-76行目
receive_force メソッドとは別に receive_message メソッドを用意します．
258-259行目
メッセージを送る方を先に説明します． Boundary の generate_force メソッドは，256行目で力を発生させた後，もし自身が床面オブジェクトだったら相手先にメッセージを送信します．メッセージは Python の辞書として作り，"type" キーに対する値を "floor_hit"，"y" キーに対する値を y 座標とします．
87-91行目
PointMass に戻ります．メッセージの処理は，pygame のイベントの処理に倣って作ります．message_list からメッセージを一つずつ取り出し，type が "floor_hit" で，かつ自身が下向きの速度成分を持っていたら，メッセージの "y" キーの値を使って自身の y 座標を更新します．
この沈み込み問題は，もっと根本的に解決できないんですか?
運動計算に拘束条件として組み込むのが正しいアプローチだと思います．

実装するのは簡単ではありませんが，基本的な考え方だけ説明しておきます．

質点は壁の内部に入り込まないという拘束条件を満たしているべきです．この拘束を満たさない場合のペナルティとなる値 (例えば侵入量の二乗) を定義し，運動計算の都度，その値を最小化するように位置の調整を行うという方法が考えられます．

同様の拘束条件は，壁との間だけではなく質点どうしの間でも満たされるべきです．(現在のプログラムはこれを全く考慮していません．複数の FixedPointMass をすり鉢状に並べて，そこに PointMass を落としてみたりするとよくわかります)．

そのほかにも，2 つの質点が剛体リンクでつながっているとか，伸び縮みしない紐で結ばれているといった状況も，質点の間の距離に関する拘束条件として表現できます．

このように複数の運動体の間に複数の拘束条件が存在する場合は，全体のペナルティが最小になるように質点系全体の位置を最適化することが必要になります．例えば非線形最小二乗問題として定式化されます．

Spring や CollisionResolver などと比べると，なんだか PointMass だけ構造が妙に複雑なように思うのですが?
はい，PointMass は抱えている責務が多すぎます．本来は複数のクラスに分割すべきです．

PointMass は力発生と運動計算という 2 つの責務を同時に担っています．物理学的に見ても質量には重力質量と慣性質量という 2 つの側面があり，本来は別物なのに一致しているものと理解されます．

少なくともこれら 2 つの責務は別のクラス (例えば GravityGenerator と MovingPoint など) に分け，PointMass クラスはそれぞれのオブジェクトへの参照を保持して (オブジェクト合成)，処理を振り分けるのが望ましい姿だと思います．

別の方針として，PointMass をこれら 2 つのクラスを継承することにより定義することも考えられます．しかし，複数のクラスを継承する多重継承(multiple inheritance) は一般に設計が難しく，下手に手を出さない方がよいとされています．Java のように多重継承を禁止している言語もあります．

7.5. 演習

例題7-1

独自の LineDrawer を定義して，Spring の表示が点線になるようにしてください．例えば 2 つの端点を結ぶ線分を適当な個数の区間に分割し，各区間につき 1 個の円を表示するなどの方法を考えるとよいです．

例題7-2

重力の無い世界で，バネ定数 \(k\) のバネの一端を固定し，もう一端につないだ質量 \(m\) の質点に，正弦波の強制振動力を加えることを考えます．強制振動の角周波数が共振角周波数 \(\sqrt{k / m}\) に近い場合と遠い場合の挙動をシミュレーションで確認してください．

(ヒント: 角周波数 \(\omega\) の単振動は，フレーム番号を \(n\), フレームあたりの時間ステップを \(\Delta t\) とすると \(\cos(\omega n \Delta t)\) で表せます．振幅は適宜設定してください)

例題7-3

CollisionResolver を継承し，累計衝突回数を画面の左上に数字で表示する機能を持つ衝突解決クラスを作成してください．

(ヒント: 例えば，累計衝突回数を保持する属性を追加し， generate_force を上書きしてそれをカウントアップするとともに， draw を上書きしてそれを表示するといった方法が考えられます)

例題7-4

Spring クラスの実装を参照しながら，2 つの質点の間に距離の二乗に反比例する引力や斥力を発生させる Actor のクラスを作成してください．