Previous | ToC | Next |
赤木: というわけで、プログラムとかの解説ね。
学生S: 随分時間がたった(作者注:作者がこれを書いているのは 2005/7 で単行本の発売直前)ので何をしたのかすっかり忘れたんですが、、、
赤木: まあ、忘れたところで新鮮な眼でみればわかりやすい解説になるんじゃない?
学生S: えー、でも、印税も入らない(実話)のにやってらんないですよ。
赤木: うるさいわね、とりあえずリスト出して。
学生S: これですね:
# octobert.rb : numerical integration for three bodies require 'vector.rb' include Math
def rotate2d(x,deg) theta=deg/180.0*(Math::PI) s = sin(theta) c = cos(theta) newx= x[0]*c-x[1]*s newy= x[0]*s+x[1]*c [newx,newy].to_v end 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 rotate(deg) @pos=rotate2d(@pos,deg) @vel=rotate2d(@vel,deg) 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 saturnsetup(saturnecc=0.0555,jupiterecc=0.0485, saturnomega=80) 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[2].convertfromelements(9.5549/5.2026,saturnecc) b.ba[2].rotate(saturnomega) b.ba[2].mass=2.8589e-4 b.ba[1].convertfromelements(1.0,jupiterecc) b.ba[1].mass=9.5479e-4 p b b.cmadjust b end
def saturnsetupb(saturnecc=0.0555,jupiterecc=0.0485, saturnomega=80) 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[2].convertfromelements(9.5549/5.2026,saturnecc) b.ba[2].rotate(saturnomega) b.ba[2].mass=0 b.ba[1].convertfromelements(1.0,jupiterecc) b.ba[1].mass=9.5479e-4 p b b.cmadjust b end def saturnsetupc(saturnecc=0.0555,jupiterecc=0.0485, saturnomega=80) 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[2].convertfromelements(9.5549/5.2026,saturnecc) b.ba[2].rotate(saturnomega) b.ba[2].mass=2.8589e-4 b.ba[1].convertfromelements(1.0,jupiterecc) b.ba[1].mass=0 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?("october.dat") endtime=50000 tpend=0 b=saturnsetup outfile=open("october.dat","w+") dt = 0.1 b.run(dt,endtime,10,tpend){ a2=b.ba[2].contactelements(b.ba[0]) a1=b.ba[1].contactelements(b.ba[0]) outfile.print b.time, " ", b.ba[1].pos[0]," ", b.ba[1].pos[1] outfile.print " ", b.ba[2].pos[0]," ", b.ba[2].pos[1] outfile.print " ",a1[0][0]," ",a1[0][1]," ",a1[1], " ",a1[2], " ",a1[3] outfile.print " ",a2[0][0]," ",a2[0][1], " ",a2[1], " ",a2[2], " ",a2[3], "\n" } outfile.close end
unless File.exist?("octoberb.dat") endtime=50000 tpend=0 b=saturnsetupb outfile=open("octoberb.dat","w+") dt = 0.15 b.run(dt,endtime,10,tpend){ a2=b.ba[2].contactelements(b.ba[0]) a1=b.ba[1].contactelements(b.ba[0]) outfile.print b.time, " ", b.ba[1].pos[0]," ", b.ba[1].pos[1] outfile.print " ", b.ba[2].pos[0]," ", b.ba[2].pos[1] outfile.print " ",a1[0][0]," ",a1[0][1]," ",a1[1], " ",a1[2], " ",a1[3] outfile.print " ",a2[0][0]," ",a2[0][1], " ",a2[1], " ",a2[2], " ",a2[3], "\n" } outfile.close end unless File.exist?("octoberc.dat") endtime=300000 tpend=0 b=saturnsetupc outfile=open("octoberc.dat","w+") dt = 0.1 b.run(dt,endtime,100,tpend){ a2=b.ba[2].contactelements(b.ba[0]) a1=b.ba[1].contactelements(b.ba[0]) outfile.print b.time, " ", b.ba[1].pos[0]," ", b.ba[1].pos[1] outfile.print " ", b.ba[2].pos[0]," ", b.ba[2].pos[1] outfile.print " ",a1[0][0]," ",a1[0][1]," ",a1[1], " ",a1[2], " ",a1[3] outfile.print " ",a2[0][0]," ",a2[0][1], " ",a2[1], " ",a2[2], " ",a2[3], "\n" } outfile.close end
unless File.exist?("octoberb2.dat") endtime=50000 tpend=0 b=saturnsetupb(0.025) outfile=open("octoberb2.dat","w+") dt = 0.15 b.run(dt,endtime,10,tpend){ a2=b.ba[2].contactelements(b.ba[0]) a1=b.ba[1].contactelements(b.ba[0]) outfile.print b.time, " ", b.ba[1].pos[0]," ", b.ba[1].pos[1] outfile.print " ", b.ba[2].pos[0]," ", b.ba[2].pos[1] outfile.print " ",a1[0][0]," ",a1[0][1]," ",a1[1], " ",a1[2], " ",a1[3] outfile.print " ",a2[0][0]," ",a2[0][1], " ",a2[1], " ",a2[2], " ",a2[3], "\n" } outfile.close end unless File.exist?("octoberb3.dat") endtime=50000 tpend=0 b=saturnsetupb(0.0555,0.0485*2) outfile=open("octoberb3.dat","w+") dt = 0.15 b.run(dt,endtime,10,tpend){ a2=b.ba[2].contactelements(b.ba[0]) a1=b.ba[1].contactelements(b.ba[0]) outfile.print b.time, " ", b.ba[1].pos[0]," ", b.ba[1].pos[1] outfile.print " ", b.ba[2].pos[0]," ", b.ba[2].pos[1] outfile.print " ",a1[0][0]," ",a1[0][1]," ",a1[1], " ",a1[2], " ",a1[3] outfile.print " ",a2[0][0]," ",a2[0][1], " ",a2[1], " ",a2[2], " ",a2[3], "\n" } outfile.close end unless File.exist?("october-fig1.ps") printer october-fig1.ps/vcps data october.dat square viewport 0.1 1.1 0.1 1.1 expand 2 limit 0 1 0 1 reloc -0.15 0.7 label e reloc 0.45 -0.15 label t expand 1.5 limit 0 50000 0.0 0.1 box BCNTSV BCNTSV color 2 xcol 1 ycol 13 conn color 3 xcol 1 ycol 8 conn color 1 pgend psfix end unless File.exist?("october-fig2.ps") printer october-fig2.ps/vcps term data october.dat square viewport 0.1 1.1 0.1 1.1 expand 2 limit 0 1 0 1 reloc -0.15 0.65 label e\\dx reloc 0.45 -0.15 label e\\dy expand 1.5 limit -0.1 0.1 -0.1 0.1 box BCNTSV BCNTSV xcol 6 ycol 7 color 3 ptype 4 0 expand 0.05 points xcol 11 ycol 12 color 2 points color 1 pgend psfix end
unless File.exist?("october-fig3a.ps") printer october-fig3a.ps/vcps data octoberb.dat square viewport 0.1 1.1 0.1 1.1 expand 2 limit 0 1 0 1 reloc -0.15 0.7 label e reloc 0.45 -0.15 label t expand 1.5 limit 0 50000 0.0 0.1 box BCNTSV BCNTSV color 1 xcol 1 ycol 13 conn color 1 pgend psfix end unless File.exist?("october-fig3b.ps") printer october-fig3b.ps/vcps term square viewport 0.1 1.1 0.1 1.1 expand 2 limit 0 1 0 1 reloc -0.15 0.65 label e\\dx reloc 0.45 -0.15 label e\\dy expand 1.5 limit -0.15 0.05 -0.1 0.1 box BCNTSV BCNTSV data octoberb.dat color 1 ptype 4 0 expand 0.05 xcol 11 ycol 12 points data octoberb2.dat color 2 xcol 11 ycol 12 points data octoberb3.dat color 3 xcol 11 ycol 12 points pgend psfix end
unless File.exist?("october-fig3c.ps") printer october-fig3c.ps/vcps data octoberc.dat square viewport 0.1 1.1 0.1 1.1 expand 2 limit 0 1 0 1 reloc -0.15 0.7 label e reloc 0.45 -0.15 label t expand 1.5 limit 0 100000 0.0 0.1 box BCNTSV BCNTSV color 1 xcol 1 ycol 8 conn color 1 pgend psfix end unless File.exist?("october-fig3d.ps") printer october-fig3d.ps/vcps term square viewport 0.1 1.1 0.1 1.1 expand 2 limit 0 1 0 1 reloc -0.15 0.65 label e\\dx reloc 0.45 -0.15 label e\\dy expand 1.5 limit -0.1 0.1 -0.1 0.1 box BCNTSV BCNTSV data octoberc.dat color 1 ptype 4 0 expand 0.05 xcol 6 ycol 7 points pgend psfix end
unless File.exist?("october-theory.dat") b=[2.07e-4,7.04e-5] k=[0.68,0.68] a=[[0, -b[0],0, k[0]*b[0]].to_v, [b[0],0, -k[0]*b[0],0 ].to_v, [0, k[1]*b[1],0, -b[1] ].to_v, [-k[1]*b[1],0, b[1],0 ].to_v] e0=rotate2d([-0.0555,0],80) e1=[-0.0485,0] e=[e0,e1].flatten.to_v
def dedt(e,a) d= a.collect{|x|x*e}.to_v end def rk2(e,a,dt) k1 = dedt(e,a)*(0.5*dt) e + dedt(e+k1,a)*dt end
t=0 dt=10 f=open("october-theory.dat", "w+"); f.print( t, " ",e.join(" "),"\n") while t < 50000 e=rk2(e,a,dt) t+= dt f.print( t, " ",e.join(" "),"\n") end f.close end unless File.exist?("october-fig4.ps") printer october-fig4.ps/vcps term data october.dat square viewport 0.1 1.1 0.1 1.1 expand 2 limit 0 1 0 1 reloc -0.15 0.65 label e\\dx reloc 0.45 -0.15 label e\\dy expand 1.5 limit -0.1 0.1 -0.1 0.1 box BCNTSV BCNTSV xcol 6 ycol 7 color 3 ptype 4 0 expand 0.05 points xcol 11 ycol 12 color 2 points color 1 lw 3 data october-theory.dat xcol 2 ycol 3 conn xcol 4 ycol 5 conn pgend psfix end
赤木: また長いわね。
学生S: たしか、8月のになんか追加したんです。
赤木: どんなふうに違うの?
学生S: えーと、憶えてません。とりあえず diff してみます。
diff -C1 august.rb october.rb *** august.rb 2004-08-22 16:18:43.000000000 +0900 --- october.rb 2005-07-27 21:20:06.000000000 +0900 *************** *** 1,3 **** # ! # august.rb : numerical integration for three bodies # --- 1,3 ---- # ! ## octobert.rb : numerical integration for three bodies # *************** *** 6,7 **** --- 6,17 ---- include Math + + def rotate2d(x,deg) + theta=deg/180.0*(Math::PI) + s = sin(theta) + c = cos(theta) + newx= x[0]*c-x[1]*s + newy= x[0]*s+x[1]*c + [newx,newy].to_v + end + class Body *************** *** 33,34 **** --- 43,50 ---- + + def rotate(deg) + @pos=rotate2d(@pos,deg) + @vel=rotate2d(@vel,deg) + end + def fromvectortocontactelements(mass,pos,vel) *************** *** 214,216 ****
! def marssetup(jupiterascale=1,marsecc=0.093,jupiterecc=0) b = Nbody.new --- 230,264 ---- ! ! def saturnsetup(saturnecc=0.0555,jupiterecc=0.0485, saturnomega=80) ! b = Nbody.new ! b.time=0 ! b.ba= [Body.new,Body.new,Body.new] (まだ一杯出るけど以下省略)
赤木: 最初のは名前が違うだけね。次のは rotate2d、ベクトルを回転 させる関数ね。これは初期条件の設定とかで使うのかしら? 学生S: えーと、、、多分そうだと思います。呼んでるのは、、、初期 条件関係だけですね。
赤木: じゃあ、クラスはいいから実際になんかしてるところにいってみて。
学生S: 最初はこれみたいですね。
unless File.exist?("october.dat") endtime=50000 tpend=0 b=saturnsetup outfile=open("october.dat","w+") dt = 0.1 b.run(dt,endtime,10,tpend){ a2=b.ba[2].contactelements(b.ba[0]) a1=b.ba[1].contactelements(b.ba[0]) outfile.print b.time, " ", b.ba[1].pos[0]," ", b.ba[1].pos[1] outfile.print " ", b.ba[2].pos[0]," ", b.ba[2].pos[1] outfile.print " ",a1[0][0]," ",a1[0][1]," ",a1[1], " ",a1[2], " ",a1[3] outfile.print " ",a2[0][0]," ",a2[0][1], " ",a2[1], " ",a2[2], " ",a2[3], "\n" } outfile.close end
えーと、何してるんでしょ? october.dat っていうファイルになんか書いてます。
赤木: 作ったファイルはどうしてるの?
学生S: それはこっちですね。
unless File.exist?("october-fig1.ps") printer october-fig1.ps/vcps data october.dat square viewport 0.1 1.1 0.1 1.1 expand 2 limit 0 1 0 1 reloc -0.15 0.7 label e reloc 0.45 -0.15 label t expand 1.5 limit 0 50000 0.0 0.1 box BCNTSV BCNTSV color 2 xcol 1 ycol 13 conn color 3 xcol 1 ycol 8 conn color 1 pgend psfix end
軸を書いて 1 列目を X 座標、 13 列目と 8 列目を Y 座標にしてグラフ書いてます。
赤木: 8月のではファイルなんか作ってなかったわよね?
学生S: えーと(と 8 月のプログラムを見る)、そうですね。これは、 計算に時間がかかるのでグラフの形とかちょっと変える度に再計算するのが嫌 になって、計算結果はファイルにしまうようにしたんだったと思います。
赤木: 普通は計算して結果書くのとグラフ書くのは別のプログラムを 使うとか、そもそもグラフ書くのは gnuplot とか使うのが普通だから、そっ ちのほうがよくあるやり方よね。
学生S: それはそうだと思います。でも、1つのプログラムで、データ ファイルとか作らないで済むならそっちのほうがわかりやすいと思います。
赤木: 私もそう思うわ。昔は MS-BASIC ではそういうのが簡単に出来 たんだけど、最近は結構そういうことに便利に使える言語とかないのよね。 MATLAB とか Mathematica とか IDL とか使えばもちろんできるんだけど、フ リーソフトじゃないし、プログラム言語としてはちょっとどうかなというとこ ろもあるし。標準的なグラフィックスがあって、で、それが例えば紙用の PS ファイルと画面に直接書くのと両方に対応してくれてるようなのってなかなか ないのよね。
まあ、それで作者は Ruby/PGPLOT を使って、さらにそれに妙な改造をしてる わけだけど。自分しか使わないプログラムってのもねえ、、、
学生S: それを使わされている僕は何か?っていう話もあります。
赤木: 君はどうせ作者の分身だからいいのよ。読者がアレよね。
学生S: でも、こんなところ読んでる読者っているんですか?
赤木: どうなのかしらね、、、と、それはともかく、後半はグラフを 書いてるってのはいいとして、前半は何してるの?
学生S: saturnsetup というのは初期条件作ってるんだと思います。中 身は
def saturnsetup(saturnecc=0.0555,jupiterecc=0.0485, saturnomega=80) 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[2].convertfromelements(9.5549/5.2026,saturnecc) b.ba[2].rotate(saturnomega) b.ba[2].mass=2.8589e-4 b.ba[1].convertfromelements(1.0,jupiterecc) b.ba[1].mass=9.5479e-4 p b b.cmadjust b end
ですから、3粒子のシステムを設定して、最初の(粒子0)は質量1で原点、最後 のは土星ということで軌道長半径と離心率から位置、速度を計算して、それか ら 80 度回してるんですね。これは本文のほうの設定通り。で、その次に木星 の位置、速度を同様に設定します。
赤木: 粒子 2 を 1 の前に設定するのはなんか深いわけがあるの?
学生S: 憶えてないですが、別にないと思います。
赤木: うーん、読むとわかりにくいわね。
学生S: そうですけど、結果は同じだし間違いではないんだからかまわ ないと思いますが。
赤木: そういってしまえばそうだけど、ほら、こんな余計な説明がい るわけじゃない。そういうのがないようにプログラムを書くのがほんとよね。
学生S: そうですか?動いて正しい答が出ればいい、っていうか、それ 以上になんかしている余裕なかったんです。
赤木: まあいいわ。 run の中身は 8 月と同じよね。ブロックを渡し て惑星とかの位置をファイルに書いてるわけね。 グラフにしてるのは時刻と、 a1[1] とか a2[1] とかね。これは結局なに?
学生S: contactelement が返してますね。これは、、、 fromvectortocontactelement が返すものそのままで、えーと、離心率です。 第一要素が離心率ベクトル、次が離心率、その次が軌道長半径で最後が近点引 数という順番で配列にいれたものを返すので。
赤木: えーと、待ってね、これ、返す配列の最初の要素はまた配列と いうかベクトルよね。後の要素はただの実数ね。
学生S: そうです。 Ruby の配列って別にそれぞれの要素の中身はなん でもいいですから。
赤木: あ、そうね。でも、それって配列じゃないわよね。
学生S: Ruby の配列って、名前が配列で番号でアクセスもできるって いうだけで実態は Lisp のリストみたいなものですから、それだけで構造体と かクラスみたいに使えちゃいます。どういう時にクラスを作ってどういう時に は配列ですますべきかっていうのは僕は書いててよくわからないです。
赤木: 指導原理は簡単よね。「楽にできるほう」計算速度とか問題な ら話は違うけど、そんなのではないでしょ?
学生S: そうですけど、どっちが楽かってやってみないとわからないじゃ ないですか。
赤木: それはそういうものよ。色々やってみるしか方法ないと思うわ。
学生S: うーん、、、そんなものですか?
赤木: と思うわ。「良いやり方」を人からおそわるのって難しいの。
学生S: 、、、、
赤木: まあ、先に進みましょ。後のグラフのプログラムは?
学生S: えーと、基本的には初期条件が違うとかプロットするものが違 うとかぐらいで殆どおんなじです。説明省略では駄目でしょうか?
赤木: うーん、書いてないのここだけじゃないしね、とりあえず他の 章のも出しておかないとね。
学生S: では、今回はここまでで。
赤木: お疲れさま。
Previous | ToC | Next |