Previous ToC Next

5. 「とんでる力学」2004年 10 月分

5.1. プログラムと解説

赤木: というわけで、プログラムとかの解説ね。

学生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