Previous | ToC | Next |
赤木: さて、問題の 11月分よね。
学生S: です。えーと、その、今回スキップ、というわけにいきません?
赤木: いかないわよね。てゆーか、なんで?
学生S: つまりですね、本文のほうみてもらえばわからないように、こ の回で使ったプログラム、作者が大昔に書いたものなんですよね、、、
赤木: えーと、まあ、そこは後でいいから君が書いたところは?
学生S: えーと、、、これですね
# plot.rb : postprocessing of three body integration 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) 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]) p e p eabs p a p argp [e,eabs,a,argp] end def contactelements(ref) fromvectortocontactelements(@mass+ref.mass, @pos-ref.pos,@vel-ref.vel) end def origincontactelements p @mass p @pos p @vel fromvectortocontactelements(1.0,@pos,@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
def xyplot @ba.each{|p| reloc p.pos[0] p.pos[1] dot } end def taplot p "taplot called" @ba.each{|p| p p a=p.origincontactelements p a reloc @time a[2] dot } end def teplot p = @ba[1] a=p.origincontactelements reloc @time a[1] dot end
def teplot2 for i in 0..(@ba.size-1) p = @ba[i] a=p.origincontactelements reloc @time a[1] color i+1 dot end color 1 end def teplot3(i) if i > 0 i=i-1 p = @ba[i] a=p.origincontactelements reloc @time a[1] color i+1 dot color 1 end end end
def vcirc(r) sqrt(1.0/r) end def plot(deltar,nbody,tend) integrator=open("|solar3","r+"); mass = 1e-5 r0=1-deltar r1=1 r2=1+deltar mass0=0 dtplot = 8 dtplot = tend/1000 if dtplot < tend/1000 integrator.puts <<-EOD 0 0 #{tend} 0.1 #{dtplot} #{nbody} #{mass} 1e-4 #{r0} 0 0 0 #{vcirc(r0)} 0 #{mass} 1e-4 #{r1} 0 0 0 #{vcirc(r1)} 0 #{mass} 1e-4 #{r2} 0 0 0 #{vcirc(r2)} 0 EOD
b = Nbody.new while s = integrator.gets sa = s.split; if sa[0] == "singlemass" nbody = sa[2].to_i b.ba = Array.new(nbody){Body.new} p b p b.ba for i in 0..(nbody-1) b.ba[i].mass = sa[3+i].to_f end end if sa[0] == "single" b.time = sa[1].to_f for i in 0..(nbody-1) for j in 0..2 b.ba[i].pos[j] = sa[2+i*6+j].to_f b.ba[i].vel[j] = sa[5+i*6+j].to_f end end b.xyplot end if sa[0] == "steps" print("t=", b.time, " ",s) end end end def plot2(deltar,nbody,tend,plotfuncstring) integrator=open("|solar3","r+"); mass = 1e-5 r0=1-deltar r1=1 r2=1+deltar mass0=0 dtplot = 32 dtplot = tend/1000 if dtplot < tend/1000 integrator.puts <<-EOD 0 0 #{tend} 0.05 #{dtplot} #{nbody} #{mass} 1e-4 #{r0} 0 0 0 #{vcirc(r0)} 0 #{mass} 1e-4 #{r1} 0 0 0 #{vcirc(r1)} 0 #{mass} 1e-4 #{r2} 0 0 0 #{vcirc(r2)} 0 EOD
b = Nbody.new while s = integrator.gets sa = s.split; if sa[0] == "singlemass" nbody = sa[2].to_i b.ba = Array.new(nbody){Body.new} p b p b.ba for i in 0..(nbody-1) b.ba[i].mass = sa[3+i].to_f end end if sa[0] == "single" b.time = sa[1].to_f for i in 0..(nbody-1) for j in 0..2 b.ba[i].pos[j] = sa[2+i*6+j].to_f b.ba[i].vel[j] = sa[5+i*6+j].to_f end b.ba[i].pos = b.ba[i].pos.to_v b.ba[i].vel = b.ba[i].vel.to_v end eval(plotfuncstring) end if sa[0] == "steps" print("t=", b.time, " ",s) end end end unless File.exist?("november-fig1.ps") printer november-fig1.ps/vcps square expand 1.0
limit -1.3 1.3 -1.3 1.3 viewport 0 0.45 0.55 1 box BCTSV BCNTSV ylabel y expand 0.5 ptype 8 0 plot(0.06,3,5000) limit -1.3 1.3 -1.3 1.3 viewport 1.0 2.0 0 1 expand 1 box BCTSV BCTSV expand 0.5 ptype 8 0 plot(0.08,3,50000)
limit -1.3 1.3 -1.3 1.3 viewport -1 0 -1 0 expand 1 box BCNTSV BCNTSV xlabel x ylabel y expand 0.5 ptype 8 0 plot(0.08,3,60000) limit -1.3 1.3 -1.3 1.3 viewport 1.0 2.0 0 1 expand 1 box BCNTSV BCTSV xlabel x expand 0.5 ptype 8 0 plot(0.1,3,900000) end
unless File.exist?("november-fig1b.ps") printer november-fig1b.ps/vcps square expand 1.0 limit -1.3 1.3 -1.3 1.3 viewport 0 0.45 0.55 1 box BCTSV BCNTSV ylabel y expand 0.5 ptype 8 0 plot(0.06,2,5000)
limit -1.3 1.3 -1.3 1.3 viewport 1.0 2.0 0 1 expand 1 box BCTSV BCTSV expand 0.5 ptype 8 0 plot(0.07,2,50000) limit -1.3 1.3 -1.3 1.3 viewport -1 0 -1 0 expand 1 box BCNTSV BCNTSV xlabel x ylabel y expand 0.5 ptype 8 0 plot(0.08,2,500000)
limit -1.3 1.3 -1.3 1.3 viewport 1.0 2.0 0 1 expand 1 box BCNTSV BCTSV xlabel x expand 0.5 ptype 8 0 plot(0.1,2,900000) end unless File.exist?("november-fig2.ps") printer november-fig2.ps/vcps square
expand 1.0 limit 0 2000 0.9 1.1 viewport 0 1 0.75 1 box BCNTSV BCNTSV ylabel a expand 0.5 ptype 8 0 plot2(0.06,3,2000,"b.taplot") expand 1.0 limit 0 60000 0.88 1.12 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a expand 0.5 ptype 8 0 plot2(0.08,3,60000,"b.taplot")
expand 1.0 limit 0 900000 0.86 1.14 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a xlabel t expand 0.5 ptype 8 0 plot2(0.1,3,900000,"b.taplot") end
unless File.exist?("november-fig2b.ps") printer november-fig2b.ps/vcps square expand 1.0 limit 0 2000 0.9 1.1 viewport 0 1 0.75 1 box BCNTSV BCNTSV ylabel a expand 0.5 ptype 8 0 plot2(0.06,2,2000,"b.taplot")
expand 1.0 limit 0 60000 0.88 1.12 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a expand 0.5 ptype 8 0 plot2(0.08,2,60000,"b.taplot") expand 1.0 limit 0 900000 0.86 1.14 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a xlabel t expand 0.5 ptype 8 0 plot2(0.1,2,900000,"b.taplot")
end unless File.exist?("november-fig2c.ps") printer november-fig2c.ps/vcps square
expand 1.0 limit 0 4000 0.9 1.1 viewport 0 1 0.8 1 box BCNTSV BCNTSV ylabel a expand 0.5 ptype 8 0 plot2(0.07,3,4000,"b.taplot") expand 1.0 limit 0 60000 0.88 1.12 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a expand 0.5 ptype 8 0 plot2(0.08,3,60000,"b.taplot")
expand 1.0 limit 0 100000 0.86 1.14 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a xlabel t expand 0.5 ptype 8 0 plot2(0.09,3,100000,"b.taplot") expand 1.0 limit 0 900000 0.86 1.14 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a xlabel t expand 0.5 ptype 8 0 plot2(0.1,3,900000,"b.taplot")
expand 1.0 limit 0 900000 0.86 1.14 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a xlabel t expand 0.5 ptype 8 0 plot2(0.11,3,1800000,"b.taplot") end
unless File.exist?("november-fig3.ps") printer november-fig3.ps/vcps square expand 1.0 viewport 0.1 1.1 0.55 1 limit 0 60000 -0.001 0.02 box BCTSV BCNTSV expand 1.5 reloc -8000 0.012 label e ptype 8 0 expand 0.5 plot2(0.08,3,60000,"b.teplot") expand 1
expand 1.0 viewport 0 1 -1 0 limit 0 60000 -0.001 0.02 box BCNTSV BCNTSV expand 1.5 reloc -8000 0.012 label e xlabel t ptype 8 0 expand 0.5 plot2(0.08,2,60000,"b.teplot") expand 1 end unless File.exist?("november-fig3b.ps") printer november-fig3b.ps/vcps square
expand 1.0 viewport 0.1 1.1 0.55 1 limit 0 60000 -0.001 0.02 box BCTSV BCNTSV expand 1.5 reloc -8000 0.012 label e ptype 8 0 expand 0.5 plot2(0.08,3,60000,"b.teplot2") expand 1 expand 1.0 viewport 0 1 -1 0 limit 0 60000 -0.001 0.02 box BCNTSV BCNTSV expand 1.5 reloc -8000 0.012 label e xlabel t ptype 8 0 expand 0.5 plot2(0.08,2,60000,"b.teplot2") expand 1 end
unless File.exist?("november-fig3c.ps") printer november-fig3c.ps/vcps term square expand 1.0 viewport 0.1 1.1 0.55 1 limit 0 900000 -0.001 0.02 box BCTSV BCNTSV expand 1.5 reloc -120000 0.012 label e ptype 8 0 expand 0.5 plot2(0.1,3,900000,"b.teplot2") expand 1
expand 1.0 viewport 0 1 -1 0 limit 0 900000 -0.001 0.02 box BCNTSV BCNTSV expand 1.5 reloc -120000 0.012 label e xlabel t ptype 8 0 expand 0.5 plot2(0.1,2,900000,"b.teplot2") expand 1 end
赤木: えーと、ちょっと待って、本文で使った図は3枚だけよね?
学生S: ですね。 november-fig1.ps と、後は 2 と 3c です。
赤木: その辺、使ってないのを消してみて。
学生S: ちょっと待って下さい。こんな感じかな?
# plot.rb : postprocessing of three body integration 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) 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]) p e p eabs p a p argp [e,eabs,a,argp] end def contactelements(ref) fromvectortocontactelements(@mass+ref.mass, @pos-ref.pos,@vel-ref.vel) end def origincontactelements p @mass p @pos p @vel fromvectortocontactelements(1.0,@pos,@vel) end
end class Nbody attr_accessor :ba, :time
def initialize @ba=[] end def plot @ba.each{|p| relocate(*p.pos) dot } end
def xyplot @ba.each{|p| reloc p.pos[0] p.pos[1] dot } end def taplot p "taplot called" @ba.each{|p| p p a=p.origincontactelements p a reloc @time a[2] dot } end
def teplot2 for i in 0..(@ba.size-1) p = @ba[i] a=p.origincontactelements reloc @time a[1] color i+1 dot end color 1 end end
def vcirc(r) sqrt(1.0/r) end def plot(deltar,nbody,tend) integrator=open("|solar3","r+"); mass = 1e-5 r0=1-deltar r1=1 r2=1+deltar mass0=0 dtplot = 8 dtplot = tend/1000 if dtplot < tend/1000 integrator.puts <<-EOD 0 0 #{tend} 0.1 #{dtplot} #{nbody} #{mass} 1e-4 #{r0} 0 0 0 #{vcirc(r0)} 0 #{mass} 1e-4 #{r1} 0 0 0 #{vcirc(r1)} 0 #{mass} 1e-4 #{r2} 0 0 0 #{vcirc(r2)} 0 EOD
b = Nbody.new while s = integrator.gets sa = s.split; if sa[0] == "singlemass" nbody = sa[2].to_i b.ba = Array.new(nbody){Body.new} p b p b.ba for i in 0..(nbody-1) b.ba[i].mass = sa[3+i].to_f end end if sa[0] == "single" b.time = sa[1].to_f for i in 0..(nbody-1) for j in 0..2 b.ba[i].pos[j] = sa[2+i*6+j].to_f b.ba[i].vel[j] = sa[5+i*6+j].to_f end end b.xyplot end if sa[0] == "steps" print("t=", b.time, " ",s) end end end def plot2(deltar,nbody,tend,plotfuncstring) integrator=open("|solar3","r+"); mass = 1e-5 r0=1-deltar r1=1 r2=1+deltar mass0=0 dtplot = 32 dtplot = tend/1000 if dtplot < tend/1000 integrator.puts <<-EOD 0 0 #{tend} 0.05 #{dtplot} #{nbody} #{mass} 1e-4 #{r0} 0 0 0 #{vcirc(r0)} 0 #{mass} 1e-4 #{r1} 0 0 0 #{vcirc(r1)} 0 #{mass} 1e-4 #{r2} 0 0 0 #{vcirc(r2)} 0 EOD
b = Nbody.new while s = integrator.gets sa = s.split; if sa[0] == "singlemass" nbody = sa[2].to_i b.ba = Array.new(nbody){Body.new} p b p b.ba for i in 0..(nbody-1) b.ba[i].mass = sa[3+i].to_f end end if sa[0] == "single" b.time = sa[1].to_f for i in 0..(nbody-1) for j in 0..2 b.ba[i].pos[j] = sa[2+i*6+j].to_f b.ba[i].vel[j] = sa[5+i*6+j].to_f end b.ba[i].pos = b.ba[i].pos.to_v b.ba[i].vel = b.ba[i].vel.to_v end eval(plotfuncstring) end if sa[0] == "steps" print("t=", b.time, " ",s) end end end unless File.exist?("november-fig1.ps") printer november-fig1.ps/vcps square expand 1.0
limit -1.3 1.3 -1.3 1.3 viewport 0 0.45 0.55 1 box BCTSV BCNTSV ylabel y expand 0.5 ptype 8 0 plot(0.06,3,5000) limit -1.3 1.3 -1.3 1.3 viewport 1.0 2.0 0 1 expand 1 box BCTSV BCTSV expand 0.5 ptype 8 0 plot(0.08,3,50000)
limit -1.3 1.3 -1.3 1.3 viewport -1 0 -1 0 expand 1 box BCNTSV BCNTSV xlabel x ylabel y expand 0.5 ptype 8 0 plot(0.08,3,60000) limit -1.3 1.3 -1.3 1.3 viewport 1.0 2.0 0 1 expand 1 box BCNTSV BCTSV xlabel x expand 0.5 ptype 8 0 plot(0.1,3,900000) pgend psfix end
unless File.exist?("november-fig2.ps") printer november-fig2.ps/vcps square expand 1.0 limit 0 2000 0.9 1.1 viewport 0 1 0.75 1 box BCNTSV BCNTSV ylabel a expand 0.5 ptype 8 0 plot2(0.06,3,2000,"b.taplot")
expand 1.0 limit 0 60000 0.88 1.12 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a expand 0.5 ptype 8 0 plot2(0.08,3,60000,"b.taplot") expand 1.0 limit 0 900000 0.86 1.14 viewport 0 1 -1.25 -0.25 box BCNTSV BCNTSV ylabel a xlabel t expand 0.5 ptype 8 0 plot2(0.1,3,900000,"b.taplot") pgend psfix
end unless File.exist?("november-fig3c.ps") printer november-fig3c.ps/vcps term square
expand 1.0 viewport 0.1 1.1 0.55 1 limit 0 900000 -0.001 0.02 box BCTSV BCNTSV expand 1.5 reloc -120000 0.012 label e ptype 8 0 expand 0.5 plot2(0.1,3,900000,"b.teplot2") expand 1 expand 1.0 viewport 0 1 -1 0 limit 0 900000 -0.001 0.02 box BCNTSV BCNTSV expand 1.5 reloc -120000 0.012 label e xlabel t ptype 8 0 expand 0.5 plot2(0.1,2,900000,"b.teplot2") expand 1 pgend psfix end
赤木: 長さ半分ね。ちょっとはいいかしら。最初の図を書いてるところは?
学生S: これですね。
unless File.exist?("november-fig1.ps") printer november-fig1.ps/vcps square expand 1.0 limit -1.3 1.3 -1.3 1.3 viewport 0 0.45 0.55 1 box BCTSV BCNTSV ylabel y expand 0.5 ptype 8 0 plot(0.06,3,5000)
limit -1.3 1.3 -1.3 1.3 viewport 1.0 2.0 0 1 expand 1 box BCTSV BCTSV expand 0.5 ptype 8 0 plot(0.08,3,50000) limit -1.3 1.3 -1.3 1.3 viewport -1 0 -1 0 expand 1 box BCNTSV BCNTSV xlabel x ylabel y expand 0.5 ptype 8 0 plot(0.08,3,60000)
limit -1.3 1.3 -1.3 1.3 viewport 1.0 2.0 0 1 expand 1 box BCNTSV BCTSV xlabel x expand 0.5 ptype 8 0 plot(0.1,3,900000) pgend psfix end
枠書いたあと plot っていうのを呼んでるだけみたいですね。
赤木: そうね。 plot の中身は?
学生S: これです:
def plot(deltar,nbody,tend) integrator=open("|solar3","r+"); mass = 1e-5 r0=1-deltar r1=1 r2=1+deltar mass0=0 dtplot = 8 dtplot = tend/1000 if dtplot < tend/1000 integrator.puts <<-EOD 0 0 #{tend} 0.1 #{dtplot} #{nbody} #{mass} 1e-4 #{r0} 0 0 0 #{vcirc(r0)} 0 #{mass} 1e-4 #{r1} 0 0 0 #{vcirc(r1)} 0 #{mass} 1e-4 #{r2} 0 0 0 #{vcirc(r2)} 0 EOD b = Nbody.new while s = integrator.gets sa = s.split; if sa[0] == "singlemass" nbody = sa[2].to_i b.ba = Array.new(nbody){Body.new} p b p b.ba for i in 0..(nbody-1) b.ba[i].mass = sa[3+i].to_f end end if sa[0] == "single" b.time = sa[1].to_f for i in 0..(nbody-1) for j in 0..2 b.ba[i].pos[j] = sa[2+i*6+j].to_f b.ba[i].vel[j] = sa[5+i*6+j].to_f end end b.xyplot end if sa[0] == "steps" print("t=", b.time, " ",s) end end end
赤木: 最初の
integrator=open("|solar3","r+");
がわからないわよね。
学生S: これはサブプロセスを作るというやつです。 solar3 という別 の、ここでは作者が書いた C のプログラムが実行されて、その入力と出力が integrator の書き込みと読み出しになるものです。 UNIX のコマンドでのパイプと同 様ですが、双方向なのがちょっと違います。
赤木: あ、だからその後の
integrator.puts <<-EOD 0 0 #{tend} 0.1 #{dtplot} #{nbody} #{mass} 1e-4 #{r0} 0 0 0 #{vcirc(r0)} 0 #{mass} 1e-4 #{r1} 0 0 0 #{vcirc(r1)} 0 #{mass} 1e-4 #{r2} 0 0 0 #{vcirc(r2)} 0 EOD
で solar3 の入力を与えてるわけね。なるほど。
学生S: そうです。その後の while では solar3 に出力を受け取って グラフを書くわけです。 solar3 の出力は
E, |DE/E0|, PE, VR: -1.50403e-05 0.000000e+00 -3.00764e-05 -0.499931 am: 0 0 2.999098985582e-05 singlemass nb= 3 1e-05 1e-05 1e-05 single 0 0.931199 0.128523 0 -0.140592 1.02172 0 0.992198 0.124675 0 -0.124674 0.992199 0 1.05303 0.121146 0 -0.111438 0.964934 0 steps = 0 E, |DE/E0|, PE, VR: -1.50403e-05 1.135362e-13 -3.00763e-05 -0.499929 single 8 -0.802115 0.485464 0 -0.531287 -0.884667 0 -0.211358 0.977229 0 -0.977314 -0.211439 0 0.511369 0.929056 0 -0.853098 0.46786 0 steps = 123 E, |DE/E0|, PE, VR: -1.50403e-05 1.485657e-13 -3.00974e-05 -0.500279 single 16 0.402918 -0.8415 0 0.936139 0.445985 0 -0.933565 -0.356606 0 0.357279 -0.934403 0 -0.511405 0.929818 0 -0.853205 -0.466523 0 steps = 251 E, |DE/E0|, PE, VR: -1.50403e-05 1.280662e-13 -3.01484e-05 -0.501125 single 24 0.137574 0.929347 0 -1.01865 0.150961 0 0.495079 -0.867864 0 0.86919 0.495659 0 -1.06268 0.0716562 0 -0.0692701 -0.966385 0
で、僕のプログラムでは singlemass で始まる行の 4, 5, 6 列目が粒子の質 量、 single で始まる行の 3 列目以降が粒子の位置、速度がその順番でとなっ てます。だから、 solar3 はそういうふうに出力するんだと思います。
赤木: 読んでから何をするかっていうと、 xyplot ね。これは?
学生S: これですね。
def xyplot @ba.each{|p| reloc p.pos[0] p.pos[1] dot } end
それぞれの粒子の xy 座標で点を打つだけです。
赤木: 他の図も同じようなものかしら?
学生S: えーと、そうですね、 plot2 ってのになってますね。これは なんか怪しい文字列を引数にしてて、さっきの xyplot の代わりにその引数で
eval(plotfuncstring)
してます。
赤木: eval って?
学生S: その文字列をプログラムとして実行するものです。 Ruby はイ ンタプリタだから、そういうこともできるわけです。
赤木: で、図2では taplot、図3では teplot2 というのを呼んでるわけね。
学生S: です。これらは軌道要素に変換して時間を横軸にしてプロット するだけです。
赤木: まあ、君が書いたところはそれでいいわ。 solar3 は?
学生S: すみません、これはちょっとまたの機会に、、、
赤木: んー、そうね、他の章先にやったほうがいいわね。
Previous | ToC | Next |