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 |