| Previous | ToC | Next |
赤木: というわけで、6月に続いてプログラムとかの解説ね。
学生S: えーと、そうですね。今回結構大変なんですよね。
august.rb : numerical integration for three bodies
require 'vector.rb'
include Math
class Body
attr_accessor :mass, :pos, :vel, :acc, :pot, :withring
def initialize
@pos=[]
@vel=[]
@withring=nil
end
def froms(s)
a=s.split.collect{|x| x.to_f}
p a
n=(a.size-1)/2
@mass = a[0]
@pos = a[1..n].to_v
@vel = a[(n+1)..(2*n)].to_v
p self
end
def convertfromelements(a,e)
@pos[0]=a*(1+e)
@pos[1]=0
@vel[0]=0
@vel[1]= sqrt(2*(1/@pos[0]-0.5/a))
@pos=@pos.to_v
@vel=@vel.to_v
end
def fromvectortocontactelements(mass,pos,vel)
raise "contactelements works only for 2D vector" if pos.size != 2
r=sqrt(pos*pos)
be = -mass/r+0.5*(vel*vel)
# be = -mass/2a, a = -mass/2be
a = -mass/(be*2)
vcirc = sqrt(mass/a)
jmax = a*vcirc
j = pos[0]*vel[1]-pos[1]*vel[0]
vbyj = vel*j
vbyh=[vbyj[1],-vbyj[0]].to_v
e = vbyh-pos*(1/r)
eabs=sqrt(e*e)
argp = atan2(-e[1],-e[0])
[e,eabs,a,argp]
end
def contactelements(ref)
fromvectortocontactelements(@mass+ref.mass,
@pos-ref.pos,@vel-ref.vel)
end
#
relation between J and e
assume a=1,j=1
ra = (1+e)
E= -1/r+v^2/2
va = sqrt(1/ra - 1/2a)
def accumulate_force(otherbody)
dr= @pos-otherbody.pos
r2 = dr*dr
r = sqrt(r2)
r3 = r2*r
mfact=@mass/r3
@acc += dr*(-otherbody.mass/r3)
otherbody.acc += dr*(@mass/r3)
@pot -= otherbody.mass/r
otherbody.pot -= @mass/r
end
def forcefromring
dummyp = Body.new
dummyp.mass=@ringmass/@ringndivide
dummyp.acc=@acc
dummyp.pot=@pot
dt = 2*Math::PI/@ringndivide
for i in 0..(@ringndivide-1)
theta = dt*i
dummyp.pos[0]=@ringradius*cos(theta)
dummyp.pos[1]=@ringradius*sin(theta)
accumulate_force(dummyp)
end
end
def ringinit(radius,mass,ndivide)
@ringradius,@ringmass,@ringndivide = radius,mass,ndivide
@withring=1
end
def makenull
null = ([0.0]*(@pos.size)).to_v
end
end
class Nbody
attr_accessor :ba, :time
def initialize
@ba=[]
end
def calculate_force
@ba.each{|x|
x.pot=0; x.acc = x.makenull
x.forcefromring if x.withring
}
for i in 0..(@ba.size-2)
for j in (i+1)..(@ba.size-1)
@ba[i].accumulate_force(@ba[j])
end
end
end
def leapfrog(first_call,dt)
if first_call
calculate_force
first_call=nil
end
@ba.each{|p| p.vel += p.acc*0.5*dt}
@ba.each{|p| p.pos += p.vel*dt}
calculate_force
@ba.each{|p| p.vel += p.acc*0.5*dt}
end
def yoshida6b(first_call,dt)
d = [0.784513610477560e0,
0.235573213359357e0,
-1.17767998417887e0,
1.31518632068391e0]
for i in 0..2 do leapfrog(first_call,dt*d[i]) end
leapfrog(first_call,dt*d[3])
for i in 0..2 do leapfrog(first_call,dt*d[2-i]) end
end
def calculate_energy
ke=pe=0
@ba.each{|p|
ke += 0.5*p.mass*(p.vel*p.vel)
pe += p.pot*p.mass
}
pe *= 0.5
[ke+pe,ke,pe]
end
def cmadjust
cmp=@ba[0].makenull
cmv=@ba[0].makenull
cmm=0
@ba.each{|p|
mass = p.mass
cmp += p.pos*mass
cmv += p.vel*mass
cmm += mass
}
cmp *= 1/cmm
cmv *= 1/cmm
@ba.each{|p|
p.pos -= cmp
p.vel -= cmv
}
end
def plot
@ba.each{|p|
relocate(*p.pos)
dot
}
end
def run(dt,tend,interval,plotendtime=0)
steps = ((tend+dt*0.5)/dt).to_i
calculate_force
e0= calculate_energy
print "Initial energy = #{e0[0]}\n"
i=0
@time=0
first_call = 1
steps.times{
yoshida6b(first_call,dt)
@time += dt
i+= 1
if i==interval or @time < plotendtime
i=0
yield
end
}
e= calculate_energy
print "Final energy = #{e[0]}, de = #{e0[0]-e[0]}\n"
end
end
the following is a test config for circular hierarchical binary
if nil
b.ba= [Body.new,Body.new]
b.ba[0].froms("0.5 1 0.0 0 1")
b.ba[1].froms("0.5 0 0 0 0")
p b
b.cmadjust
b.ba += [Body.new]
#
# M=2, R=5, a= 2/25, v^2/=a v=sqrt(2/5)
# pe = 2/5 ke=1/5
b.ba[2].froms("1 5 0 0 0")
b.ba[2].vel[1]=sqrt(2.0/5.0)
b.cmadjust
end
def marssetup(jupiterascale=1,marsecc=0.093,jupiterecc=0)
b = Nbody.new
b.time=0
b.ba= [Body.new,Body.new,Body.new]
b.ba[0].froms("1 0.0 0.0 0 0")
b.ba[1].convertfromelements(5.203/1.524*jupiterascale,jupiterecc)
b.ba[1].mass=0.001
b.ba[2].convertfromelements(1,marsecc)
b.ba[2].mass=0.00000
p b
b.cmadjust
b
end
class Nbody
def connectorplot(dt,tconnect)
if @time <= dt or @time > tconnect
relocate(*ba[2].pos)
dot
else
draw(*ba[2].pos)
end
end
end
unless File.exist?("august-fig1b.ps")
printer august-fig1b.ps/vcps
term
endtime=10000
endtime=100
tpend=6.38
b=marssetup
square
expand 1.5
viewport 0 0.45 0.55 1
limit -1.2 1.2 -1.2 1.2
box BCTSV BCNTSV
ylabel y
expand 0.1
dt = 0.07
b.run(dt,endtime,50,tpend){
b.connectorplot(dt,tpend)
}
viewport 1.2 2.2 0 1
b=marssetup
b.ba[1].mass*10
expand 1.5
limit -1.2 1.2 -1.2 1.2
box BCTSV BCTSV
expand 0.1
b.run(dt,endtime,50,tpend){
b.connectorplot(dt,tpend)
}
viewport -1.2 -0.2 -1.2 -0.2
b=marssetup(0.5)
expand 1.5
limit -1.2 1.2 -1.2 1.2
box BCNTSV BCNTSV
xlabel x
ylabel y
expand 0.1
b.run(dt,endtime,50,tpend){
b.connectorplot(dt,tpend)
}
expand 1
viewport 1.2 2.2 0 1
b=marssetup(0.5)
b.ba[1].mass*10
expand 1.5
limit -1.2 1.2 -1.2 1.2
box BCNTSV BCTSV
xlabel x
expand 0.1
b.run(dt,endtime,50,tpend){
b.connectorplot(dt,tpend)
}
expand 1
pgend
psfix
end
unless File.exist?("august-fig2.ps")
printer august-fig2.ps/vcps
b = marssetup
square
viewport 0.1 1.1 0.1 1.1
limit 0 500 -0.005 0.015
expand 1.2
box BCNSTV BCNSTV
expand 2
xlabel T
reloc -100 0.007
label \\gh
a=b.ba[2].contactelements(b.ba[0])
p b.time
p a[3]
reloc b.time a[3]
b.run(0.07,500,5){
a=b.ba[2].contactelements(b.ba[0])
draw(b.time,a[3])
}
b = marssetup
b.ba[1].mass*=0.5
a=b.ba[2].contactelements(b.ba[0])
reloc b.time a[3]
b.run(0.07,500,5){
a=b.ba[2].contactelements(b.ba[0])
draw(b.time,a[3])
}
b = marssetup
b.ba[1].mass*=0
a=b.ba[2].contactelements(b.ba[0])
reloc b.time a[3]
b.run(0.07,500,5){
a=b.ba[2].contactelements(b.ba[0])
draw(b.time,a[3])
}
expand 1
pgend
psfix
end
unless File.exist?("august-fig3.ps")
printer august-fig3.ps/vcps
$plotlylimit = -6
$ymin = exp($plotlylimit*log(10.0))
$lymax = $plotlylimit
def logp(x,y)
p x
p y
if y > 0
ly = log10(y)
if ly > $lymax
draw x ly
$lymax = ly
end
end
end
square
viewport 0.1 1.1 0.1 1.1
limit 0 500 $plotlylimit -0.5
expand 1.2
box BCNSTV BCNSTVL
expand 2
xlabel T
ylabel \\gh
b = marssetup
a=b.ba[2].contactelements(b.ba[0])
reloc b.time $plotlylimit
b.run(0.07,500,5){
a=b.ba[2].contactelements(b.ba[0])
logp(b.time, a[3])
}
$lymax = $plotlylimit
b = marssetup
b=marssetup(0.5)
a=b.ba[2].contactelements(b.ba[0])
reloc b.time $plotlylimit
b.run(0.07,500,5){
a=b.ba[2].contactelements(b.ba[0])
logp(b.time, a[3])
}
$lymax = $plotlylimit
b = marssetup
b=marssetup(2)
a=b.ba[2].contactelements(b.ba[0])
reloc b.time $plotlylimit
b.run(0.07,500,5){
a=b.ba[2].contactelements(b.ba[0])
logp(b.time, a[3])
}
$lymax = $plotlylimit
b = marssetup
b=marssetup(4)
a=b.ba[2].contactelements(b.ba[0])
reloc b.time $plotlylimit
b.run(0.05,500,5){
a=b.ba[2].contactelements(b.ba[0])
logp(b.time, a[3])
}
end
unless File.exist?("august-fig4.ps")
printer august-fig4.ps/vcps
square
viewport 0.1 1.1 0.1 1.1
limit 0 500 -0.005 0.01
expand 1.2
box BCNSTV BCNSTV
expand 2
xlabel T
ylabel \\gh
b = marssetup(1,0.1)
a=b.ba[2].contactelements(b.ba[0])
reloc b.time 0
b.run(0.07,500,5){
a=b.ba[2].contactelements(b.ba[0])
draw b.time a[3]
}
b = marssetup(1,0.4)
a=b.ba[2].contactelements(b.ba[0])
reloc b.time 0
b.run(0.03,500,5){
a=b.ba[2].contactelements(b.ba[0])
draw b.time a[3]
}
end
unless File.exist?("august-fig5.ps")
printer august-fig5.ps/vcps
term
square
viewport 0.1 1.1 0.1 1.1
tend = 500
limit 0 tend -0.005 0.01
expand 1.2
box BCNSTV BCNSTV
expand 2
xlabel T
ylabel \\gh
color 1
b = marssetup
reloc b.time 0
b.run(0.07,tend,5){
a=b.ba[2].contactelements(b.ba[0])
draw b.time a[3]
}
b = marssetup
reloc b.time 0
b.ba[1].mass=0
b.ba[2].ringinit(5.203/1.524,0.001,32)
color 2
b.run(0.07,tend,5){
a=b.ba[2].contactelements(b.ba[0])
draw b.time a[3]
}
color 1
end
赤木: これちょっと長いどころの騒ぎじゃないわね。 500行近いわよ。
学生S: これ1つで8月号の図全部書くんだから、しょうがないと思いま す。頭から見ていってもいいんですが、どう使われるかをみないとわけがわか らないところがおおいと思いますから、実際に図を書くところから見ていき ませんか？
赤木: そうね。じゃあ、最初の図のところね。
学生S: はい、メインプログラムというか、トップレベルの、関数とか の中でなくて図1を書いてるのはここです。
unless File.exist?("august-fig1b.ps")
printer august-fig1b.ps/vcps
term
endtime=10000
endtime=100
tpend=6.38
b=marssetup
square
expand 1.5
viewport 0 0.45 0.55 1
limit -1.2 1.2 -1.2 1.2
box BCTSV BCNTSV
ylabel y
expand 0.1
dt = 0.07
b.run(dt,endtime,50,tpend){
b.connectorplot(dt,tpend)
}
viewport 1.2 2.2 0 1
b=marssetup
b.ba[1].mass*10
expand 1.5
limit -1.2 1.2 -1.2 1.2
box BCTSV BCTSV
expand 0.1
b.run(dt,endtime,50,tpend){
b.connectorplot(dt,tpend)
}
viewport -1.2 -0.2 -1.2 -0.2
b=marssetup(0.5)
expand 1.5
limit -1.2 1.2 -1.2 1.2
box BCNTSV BCNTSV
xlabel x
ylabel y
expand 0.1
b.run(dt,endtime,50,tpend){
b.connectorplot(dt,tpend)
}
expand 1
viewport 1.2 2.2 0 1
b=marssetup(0.5)
b.ba[1].mass*10
expand 1.5
limit -1.2 1.2 -1.2 1.2
box BCNTSV BCTSV
xlabel x
expand 0.1
b.run(dt,endtime,50,tpend){
b.connectorplot(dt,tpend)
}
expand 1
pgend
psfix
end
赤木: で、これが何をしているの？
学生S: 要するに、"august-fig1b.ps" という ファイルに図1を書く、というものですね。
赤木: 1b なのはなんか深いわけがあるの？
学生S: えーと、すみません、だいぶ前の話なんで忘れました。という か、これ編集から注文が来て書き直したんで変わってるんだったと思います。
赤木: プログラムのほうが残してあれば、図とかのファイルの名前は 変えなくてもいいような気がするけど、まあいいわ。最初の
unless File.exist?("august-fig1b.ps")
っていうのは、図があれば新しくは作らないというのね？
学生S: そうです。トップレベルで
unless File.exist?("foo.ps")
...
end
というのが一杯あって、それぞれ別の図を書くというふうになってます。図毎 にプログラムわけてもいいんですが、そうするとかえって複雑になるような気 がして。
赤木: まあ、原稿一つの図が一つのプログラムで全部かけるというの は便利は便利かもね。 次の
printer august-fig1b.ps/vcps
は Rongo のコマンドになってるわけね？普通の Ruby/PGPLOT ではなくて。
学生S: です。
赤木: その後いくつか変数を初期化して、、、あ、
b=marssetup
ってのは何？
学生S: これは時間積分の初期設定です。marssetup の実体は
def marssetup(jupiterascale=1,marsecc=0.093,jupiterecc=0)
b = Nbody.new
b.time=0
b.ba= [Body.new,Body.new,Body.new]
b.ba[0].froms("1 0.0 0.0 0 0")
b.ba[1].convertfromelements(5.203/1.524*jupiterascale,jupiterecc)
b.ba[1].mass=0.001
b.ba[2].convertfromelements(1,marsecc)
b.ba[2].mass=0.00000
p b
b.cmadjust
b
end
で、Nbody というクラスの変数 b を作って、時刻 @time を 0 にして粒子3個 の質量、位置、速度を初期設定します。粒子は @ba という配列にはいってい て、まず Body.new で初期化してから、粒子0(太陽)は質量 1、位置、速度は 0 で、 木星、火星(1、2)は軌道長半径と離心率から位置、速度を計算します。計算す るのは
def convertfromelements(a,e)
@pos[0]=a*(1+e)
@pos[1]=0
@vel[0]=0
@vel[1]= sqrt(2*(1/@pos[0]-0.5/a))
@pos=@pos.to_v
@vel=@vel.to_v
end
で、最初は x 軸上で、遠日点であるとして位置、速度を計算するようになっ てます。
赤木: その後しばらくは画面の設定ね。それはいいとして
b.run(dt,endtime,50,tpend){
b.connectorplot(dt,tpend)
}
だけど、これで {} の中は？なんかいかにもこれで絵を書いてるって感じだけ ど、、、
学生S: そうですね、これは Ruby のちょっと便利な機能を使ってみま した。 Nbody#run (Nbody というクラスの run というメソッドという意味) は
def run(dt,tend,interval,plotendtime=0)
steps = ((tend+dt*0.5)/dt).to_i
calculate_force
e0= calculate_energy
print "Initial energy = #{e0[0]}\n"
i=0
@time=0
first_call = 1
steps.times{
yoshida6b(first_call,dt)
@time += dt
i+= 1
if i==interval or @time < plotendtime
i=0
yield
end
}
e= calculate_energy
print "Final energy = #{e[0]}, de = #{e0[0]-e[0]}\n"
end
というものなんですが、ここで
yield
というのがありますよね？プログラムのここにくると、 {} の中身のコードが 実行されるんです。Ruby ではこういうのをブロックとかコードブロックといいます。
赤木: えーと、関数をパラメータで渡すみたいなもの？
学生S: まあ、この場合はそうともいえますね。中でメソッド一つ呼ん でるだけですし。でも、ここでしか使わないものだといちいち関数程度しなく ても使えるから簡単でわかりやすいし。あと、これがちょっと面白いのは、 ブロックが、呼ばれた関数、この場合 run ですね、のコンテキストではなく て呼んだほう、この場合はトップレベルですが、のほうのコンテキストで実行 されるということです。つまり、呼んでほうのローカル変数とかは見えるけど 呼ばれたほうのは見えないんです。
赤木: それって不便っていうか、やりたいことと違うんではないの？
学生S: 実はこの場合にはちょっとそういうところもあるんですけど。 一般的にはこっちのほうが便利だと思います。Ruby のいろんな機能自体がこ れを使って実現されてますから。例えば each とかって使ったと思いますけど、 これは実はメソッドなわけでブロックが呼んだ側のコンテキストで実行されな いと each に渡すブロックに何も書けなくなってしまいます。
赤木: それはそうね。
学生S: で、先に、、、といっても、 run を呼んだ後はもう絵を書く 処理を終えてるだけですね。 run のほうの中身にいってみると、
まず
steps = ((tend+dt*0.5)/dt).to_i
は繰り返しの回数計算ですね。で、ここではいわゆる吉田6次公式というシン プレクティック積分公式を使うのですが、そのためには最初に加速度を計算し ておく必要があるので
calculate_force
で加速度を計算してます。これは後で説明します。次はチェック用ですが、
e0= calculate_energy
print "Initial energy = #{e0[0]}\n"
で系の全エネルギーを計算してます、後は繰り返しです。
i=0
@time=0
first_call = 1
steps.times{
yoshida6b(first_call,dt)
これが時間積分ですね。
@time += dt
i+= 1
if i==interval or @time < plotendtime
ここで、条件がややこしいのは、編集から最初の1周期は線をひいて欲しいと いう要請があったからです。それ以降は何回かに一回だけ点を打つわけです。 その判断は、問題の
yield
で呼ばれた
b.connectorplot(dt,tpend)
のほうがしているわけです。
赤木: なるほど。まあ、プログラム自体は割合素直ね。とすると後は
calculate_force
calculate_energy
yoshida6b(first_call,dt)
connectorplot(dt,tpend)
ね、順番にいってみて。
学生S: まず calculate_force ですね。
def calculate_force
@ba.each{|x|
x.pot=0; x.acc = x.makenull
x.forcefromring if x.withring
}
for i in 0..(@ba.size-2)
for j in (i+1)..(@ba.size-1)
@ba[i].accumulate_force(@ba[j])
end
end
end
うーん、見ての通りですが、、、
赤木: 最初は全粒子についてポテンシャル、加速度を 0 にしてるのね。 makenull ってのは？
学生S: 0 ベクトルを作るだけの関数です。その次に、後で木星をリン グで置き換えるのがあるのでその時にはリングからの力を計算します。これは こんなの
def forcefromring
dummyp = Body.new
dummyp.mass=@ringmass/@ringndivide
dummyp.acc=@acc
dummyp.pot=@pot
dt = 2*Math::PI/@ringndivide
for i in 0..(@ringndivide-1)
theta = dt*i
dummyp.pos[0]=@ringradius*cos(theta)
dummyp.pos[1]=@ringradius*sin(theta)
accumulate_force(dummyp)
end
end
で、円環上にダミー粒子を並べて、そこからの力を計算するというのです。
赤木: これ、呼ばれる度にダミー粒子作り直すの？
学生S: そうです。無駄は無駄なんですが、そんなに時間かからなかっ たというか、書き直すのにかかりそうな時間より速く計算終わったので。
赤木: まあ、いいわ。 accumulate_force が実際に加速度を計算して るわけね。
学生S: そうです。
def accumulate_force(otherbody)
dr= @pos-otherbody.pos
r2 = dr*dr
r = sqrt(r2)
r3 = r2*r
mfact=@mass/r3
@acc += dr*(-otherbody.mass/r3)
otherbody.acc += dr*(@mass/r3)
@pot -= otherbody.mass/r
otherbody.pot -= @mass/r
end
です。これはもう式の通りです。引数に相手の粒子をとって、そちらにも積算 します。ですから、 calculate_force では内側のループは i+1 から回すよう になってます。
calculate_energy は、
def calculate_energy
ke=pe=0
@ba.each{|p|
ke += 0.5*p.mass*(p.vel*p.vel)
pe += p.pot*p.mass
}
pe *= 0.5
[ke+pe,ke,pe]
end
これも見ての通りですが、、、あ、最後の
[ke+pe,ke,pe]
は、これらを要素にする配列を作っていて、このメソッドの返り値がこの配列 になります。こういうふうに気軽に配列とかを返せるのは便利だと思います。時間積分は
def yoshida6b(first_call,dt)
d = [0.784513610477560e0,
0.235573213359357e0,
-1.17767998417887e0,
1.31518632068391e0]
for i in 0..2 do leapfrog(first_call,dt*d[i]) end
leapfrog(first_call,dt*d[3])
for i in 0..2 do leapfrog(first_call,dt*d[2-i]) end
end
で、これは都合7回、時間刻みをいろいろ変えた leapfrog を呼ぶというもの ですね。leapfrog は
def leapfrog(first_call,dt)
if first_call
calculate_force
first_call=nil
end
@ba.each{|p| p.vel += p.acc*0.5*dt}
@ba.each{|p| p.pos += p.vel*dt}
calculate_force
@ba.each{|p| p.vel += p.acc*0.5*dt}
end
で、公式通りです。グラフを書く connectorplot は
def connectorplot(dt,tconnect)
if @time <= dt or @time > tconnect
relocate(*ba[2].pos)
dot
else
draw(*ba[2].pos)
end
end
で、指定した時刻の前なら線を引いて、後なら点を打つというだけです。
relocate(*ba[2].pos)
で、
*ba[2].pos
というのは見なれないと思いますが、これは、ここで配列が、なんというのか、 展開されて、要素が2つなので
ba[2].pos[0], ba[2].pos[1]
と書くのと同じことになります。
赤木: んー、なくても困らないような機能ね。
学生S: 困ることはないかもしれないですけど、短く書けるのはいいと 思います。
赤木: まあ、それはそうね。他の図は？
学生S: 全部 run でブロックを渡して書くというのは今回同じなので、 ここまででいいことにしませんか？
赤木: ちょっと手抜きな気がするけど、、、
学生S: でも、ほら、これ8月号分で、9月号が今週でちゃいますよね？ そっち書いて時間があったらこっちももうちょっと追加ということで、、、
赤木: そのほうがいいのは確かね。じゃあ、そうさせてもらいましょ。 あ、ルンゲ・レンツベクトルの話は？
学生S: え、でも、吉田さんの教科書とか見てもらったほうがいいと思 いますけど、、、
赤木: まね。リンクは？
学生S:これですね。
| Previous | ToC | Next |