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