🌏
万有引力のシミュレーション
地球は別に円運動をしたいわけではなく太陽にひっぱられている
ニュートンの万有引力の法則を確認する
F = 万有引力定数 * 質量1 * 質量2 * / 距離**2
万有引力定数は 6.67259 と宇宙全体で決まっているのだけど、ここでは 1 として無視する。そして太陽と地球の質量をそれぞれ 100, 2 として、
太陽 = 100.0
地球 = 2.0
距離だけを変更してみる。
f = 太陽 * 地球 / 10000**2 # => 2.0e-06
f = 太陽 * 地球 / 1000**2 # => 0.0002
f = 太陽 * 地球 / 100**2 # => 0.02
f = 太陽 * 地球 / 10**2 # => 2.0
f = 太陽 * 地球 / 1**2 # => 200.0
すると近づくにつれて引力が急激に強くなっていくのがわかる。つまり近いとめちゃくちゃ強い。遠いとめちゃくちゃ弱い。
引力がそのまま加速度になるわけではない
ここが間違えやすい。仮に太陽と地球の間で 100 の引力が生まれたとき、
太陽の加速度 = 100
地球の加速度 = 100
とするのは、地球より遥かに大きい太陽が地球と同じ速度で寄ってくることになるのでおかしい。なのでここは自分の質量で割るのが正しい。[1]
太陽の加速度 = 100 / 太陽 # => 1.0
地球の加速度 = 100 / 地球 # => 50.0
地球視点の実装手順
万有引力定数
gravitation = 0.5
太陽
attractor_location = V[400.0, 400.0]
attractor_mass = 100.0
地球
location = V[400.0, 200.0]
velocity = V[8.0, 0.0]
mass = 2.0
引力を求める(両方の視点で共通)
diff = attractor_location - location # => (0.0, 200.0)
distance = diff.length # => 200.0
distance = distance.clamp(5, 15) # => 15
strength = gravitation * mass * attractor_mass / distance**2 # => 0.4444444444444444
force = strength * diff.normalize # => (0.0, 0.4444444444444444)
地球の加速度に適用する
acceleration = force / mass # => (0.0, 0.2222222222222222)
地球を移動する
velocity += acceleration # => (8.0, 0.2222222222222222)
location += velocity # => (408.0, 200.22222222222223)
距離を制限する理由
distance = distance.clamp(5, 15) # => 15
これを入れないと地球が吹っ飛ぶから。
たとえば太陽に近くなると、ものすごい引力が生まれて太陽につっこんでしまうので、近すぎたときはそんなに近くなってないことにする。
また太陽からかなり離れてしまうと、ほとんど引力の影響を受けなくなっておもんないので、遠すぎたときはそんなに遠くなってないことにする。
万有引力定数の用途
惑星のそれぞれの質量を「音量」としたとき万有引力定数は「マスター音量」と考えるとわかりやすい。つまり 0 にすると宇宙全体の引力が消滅し、1 にすると各引力がそのまま適用される。
マスター音量のようなパラメータであれば常に 1 とすれば別になくてもよいのだけど、もし全体の引力を反発力に切り替えたい場合、符号を反転させるだけでよいという使い道もある。
相互引力の実験
中央に太陽があるわけではない
同じ質量を持った惑星が散らばっているとした場合、自然とクラスタを形成する。上のデモは画面の中央に何かがあったわけではなく、各惑星が自分以外のすべての惑星に少しずつ引っ張られた結果、なんとなく全体の中心に引き寄せられることになり、それがたまたま画面の中央だっただけ。
引力を反転するとどうなる？
同じ質量を持った惑星が反発し合うとクラスタの中心から離れていく。これも同様に、すべての惑星から少しずつ反発しただけだが、クラスタの中心を基点に反発したように感じる。
実験用の最小コード
class App < Base
def initialize
super
reset
end
def reset
@points = [window_wh.center]
@location = window_wh.center + V.from_angle(270.deg_to_rad).clamp_length_min(220)
@velocity = V.from_angle(0.deg_to_rad).clamp_length_min(8)
end
def button_down(id)
super
case id
when Gosu::KB_R
reset
end
end
def draw
super
gravitation = 0.5
# 太陽
attractor_location = @points.first
attractor_mass = 100.0
# 地球
mass = 1.0
radius = 50
diff = attractor_location - @location
distance = diff.length
distance = distance.clamp(5, 15)
strength = gravitation * mass * attractor_mass / distance**2
force = strength * diff.normalize_or_zero
@velocity += force / mass
@location += @velocity
vputs "万有引力定数: #{gravitation}"
vputs "太陽の質量: #{attractor_mass}"
vputs "地球の質量: #{mass}"
vputs "距離: #{diff.magnitude.round}"
vputs "距離(補正後): #{distance.round}"
vputs "地球の加速度: #{strength.round(2)}"
circle_draw(@location, radius, sides: 32)
line_draw(@location, @location + @velocity.clamp_length_min(radius))
circle_draw(attractor_location, 100, sides: 32)
end
def window_size_default
V[800, 800]
end
show
end
パーティクル
class Planet
attr_accessor :location
attr_accessor :velocity
attr_accessor :radius
attr_accessor :mass
def initialize
@location = $app.window_wh.center + V.rand_norm * $app.window_wh.center
@velocity = V.rand_norm
@acceleration = V.zero
@radius = 16
@mass = 1
end
def attract_to(other)
diff = other.location - @location
distance = diff.length.clamp($app.min, $app.max)
strength = $app.gravitation * @mass * other.mass / distance**2
force = strength * diff.normalize
if $app.reverse_p
force *= -1
end
force
end
def update
if $app.attractor_p
$app.attractors.each do |attractor|
@acceleration += attract_to(attractor) / @mass
end
else
$app.planets.each do |other|
next if self == other
@acceleration += attract_to(other) / @mass
end
end
@velocity += @acceleration
@velocity *= $app.damping
@velocity = @velocity.clamp_length_min(2.0) if false
@location += @velocity
@acceleration = V.zero
if $app.frame_p
if @location.x < 0 || $app.window_wh.x <= @location.x || @location.y < 0 || $app.window_wh.y <= @location.y
@velocity = -@velocity
end
end
end
def draw
$app.point_draw(@location)
end
def acceleration_add(force)
@acceleration += force / @mass
end
end
class Attractor
attr_accessor :mass
attr_accessor :radius
attr_accessor :location
attr_accessor :index
def initialize(index)
@index = index
@mass = 50
@radius = 50.0
@location = $app.window_wh.center
end
def draw
if $app.attractor_p
$app.circle_draw(location, @radius, sides: 32)
end
end
end
class App < Base
attr_accessor :attractors
attr_accessor :planets
attr_accessor :attractor_p
attr_accessor :frame_p
attr_accessor :damping
attr_accessor :gravitation
attr_accessor :reverse_p
attr_accessor :min
attr_accessor :max
def initialize
super
preset1
end
def preset1
@attractor_p = true
@frame_p = false
@damping = 0.998
@gravitation = 0.5
@reverse_p = false
@min = 5
@max = 30
@attractors = 1.times.collect.with_index { |i| Attractor.new(i) }
@points = @attractors.collect(&:location)
@planets = 100.times.collect { Planet.new }
end
def button_down(id)
super
case id
when Gosu::KB_1
preset1
when Gosu::KB_A
@attractor_p ^= true
when Gosu::KB_F
@frame_p ^= true
when Gosu::KB_R
@reverse_p ^= true
when Gosu::KB_J, Gosu::KB_K
@gravitation += (id == Gosu::KB_K ? -1 : 1) * 0.01
when Gosu::KB_H, Gosu::KB_L
@damping += (id == Gosu::KB_H ? -1 : 1) * 0.001
when Gosu::KB_UP, Gosu::KB_DOWN
@min += (id == Gosu::KB_UP ? -1 : 1) * 1
@min = [@min, @max].min
when Gosu::KB_LEFT, Gosu::KB_RIGHT
@max += (id == Gosu::KB_LEFT ? -1 : 1) * 1
@max = [@min, @max].max
end
end
def draw
super
@planets.each(&:update)
@planets.each(&:draw)
@attractors.each(&:draw)
vputs "[a] アトラクタ: #{@attractor_p}"
vputs "[r] 反転: #{@reverse_p}"
vputs "[f] 外枠: #{@frame_p}"
vputs "[h][l] 減退率: #{@damping.round(3)}"
vputs "[j][k] 引力: #{@gravitation.round(2)}"
vputs "[↑][↓] 距離最小: #{@min}"
vputs "[←][→] 距離最大: #{@max}"
end
def window_size_default
V[800, 800]
end
show
end
参照
- https://natureofcode.com/book/chapter-2-forces/#:~:text=2.9 Gravitational Attraction
- https://ja.wikipedia.org/wiki/万有引力
- https://github.com/mark-gerarts/nature-of-code/blob/master/02. Forces/Exercise 2.9 - Custom force/sketch.lisp
-
Nature of Code でもそうなっている ↩︎
Discussion