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