13. プラマーモデル
赤木 ：前回はとりあえず粒子データ読み書きだけだったから、今日は なんかそれらしい初期モデルを作りましょう。
学生C：作りましょうっていわれても、、、
赤木 ：分布関数とか簡単に書けるのはプラマーモデル Plummer model という やつね。これは、球対称で、星の密度部分が、半径とか密度の値を適当にスケー ルすると
| Previous | ToC | Next |
赤木 ：前回はとりあえず粒子データ読み書きだけだったから、今日は
なんかそれらしい初期モデルを作りましょう。
学生C：作りましょうっていわれても、、、
赤木 ：分布関数とか簡単に書けるのはプラマーモデル Plummer model という
やつね。これは、球対称で、星の密度部分が、半径とか密度の値を適当にスケー
ルすると
13. プラマーモデル
で与えられて、重力ポテンシャルが
になるのね。で、初期条件の作り方がこの
論文
の Appendix に詳しく書いてあるから、みんなまずはこれ使うの。
あとはもちろん、 Hernquist モデルとかその一般化とか NFWプロファイルと かその一般化とかあるんだけど、まあとりあえずは。
学生C：えっと、その論文読んで作るんですか？
赤木 ：あ、それでもいいけど、 Ruby のがあるから、それを Crystal で動く ように直すのでもいいわ。 Ruby のはこれ。
1:#!/usr/local/bin/ruby -w
2:
3:require "acs.rb"
4:require "acs"
5:options_text= <<-END
6:
7: Description: Plummer's Model Builder
8: Long description:
9: This program creates an N-body realization of Plummer's Model.
10: (c) 2004, Piet Hut and Jun Makino; see ACS at www.artcompsi.org
11:
12: The algorithm used is described in Aarseth, S., Henon, M., & Wielen, R.,
13: Astron. Astroph. 37, 183 (1974).
14:
15:
16: Short name:-n
17: Long name: --n_particles
18: Value type: int
19: Default value: 1
20: Variable name: n
21: Print name: N
22: Description: Number of particles
23: Long description:
24: Number of particles in a realization of Plummer's Model.
25:
26: Each particles is drawn at random from the Plummer distribution,
27: and therefore there are no correlations between the particles.
28:
29: Standard Units are used in which G = M = 1 and E = -1/4, where
30: G is the gravitational constant
31: M is the total mass of the N-body system
32: E is the total energy of the N-body system
33:
34:
35: Short name: -s
36: Long name: --seed
37: Value type: int
38: Default value: 0
39: Description: pseudorandom number seed given
40: Print name:
41: Variable name: seed
42: Long description:
43: Seed for the pseudorandom number generator. If a seed is given with
44: value zero, a preudorandom number is chosen as the value of the seed.
45: The seed value used is echoed separately from the seed value given,
46: to allow the possibility to repeat the creation of an N-body realization.
47:
48: Example:
49:
50: |gravity> kali mkplummer1.rb -n 42 -s 0
51: . . .
52: pseudorandom number seed given: 0
53: actual seed used: 1087616341
54: . . .
55: |gravity> kali mkplummer1.rb -n 42 -s 1087616341
56: . . .
57: pseudorandom number seed given: 1087616341
58: actual seed used: 1087616341
59: . . .
60:
61:
62: END
63:
64:class Body
65:
66: attr_accessor :body_id, :mass, :pos, :vel
67:
68: def ekin # kinetic energy
69: 0.5*@mass*(@vel*@vel)
70: end
71:
72: def epot(body_array) # potential energy
73: p = 0
74: body_array.each do |b|
75: unless b == self
76: r = b.pos - @pos
77: p += -@mass*b.mass/sqrt(r*r)
78: end
79: end
80: p
81: end
82:
83:end
84:
85:class NBody
86:
87: attr_accessor :time, :body
88:
89: def initialize
90: @body = []
91: end
92:
93: def ekin # kinetic energy
94: e = 0
95: @body.each{|b| e += b.ekin}
96: e
97: end
98:
99: def epot # potential energy
100: e = 0
101: @body.each{|b| e += b.epot(@body)}
102: e/2 # pairwise potentials were counted twice
103: end
104:
105: def adjust_center_of_mass
106: vel_com = pos_com = @body[0].pos*0 # null vectors of the correct length
107: @body.each do |b|
108: pos_com += b.pos*b.mass
109: vel_com += b.vel*b.mass
110: end
111: @body.each do |b|
112: b.pos -= pos_com
113: b.vel -= vel_com
114: end
115: end
116:
117: def adjust_units
118: alpha = -epot / 0.5
119: beta = ekin / 0.25
120: @body.each do |b|
121: b.pos *= alpha
122: b.vel /= sqrt(beta)
123: end
124: end
125:
126:end
127:
128:def frand(low, high)
129: low + rand * (high - low)
130:end
131:
132:def spherical(r)
133: vector = Vector.new
134: theta = acos(frand(-1, 1))
135: phi = frand(0, 2*PI)
136: vector[0] = r * sin( theta ) * cos( phi )
137: vector[1] = r * sin( theta ) * sin( phi )
138: vector[2] = r * cos( theta )
139: vector
140:end
141:
142:def mkplummer(c)
143: if c.seed == 0
144: srand
145: else
146: srand c.seed
147: end
148: nb = NBody.new
149: c.n.times do |i|
150: b = plummer_sample
151: b.mass = 1.0/c.n
152: b.body_id = i
153: nb.body.push(b)
154: end
155: nb.adjust_center_of_mass if c.n > 0
156: nb.adjust_units if c.n > 1
157: nb.acs_log(1, " actual seed used\t: #{srand}\n")
158: nb.acs_write($stdout, false, c.precision, c.add_indent)
159:end
160:
161:def plummer_sample
162: b = Body.new
163: scalefactor = 16.0 / (3.0 * PI)
164: radius = 1.0 / sqrt( rand ** (-2.0/3.0) - 1.0)
165: b.pos = spherical(radius) / scalefactor
166: x = 0.0
167: y = 0.1
168: while y > x*x*(1.0-x*x)**3.5
169: x = frand(0,1)
170: y = frand(0,0.1)
171: end
172: velocity = x * sqrt(2.0) * ( 1.0 + radius*radius)**(-0.25)
173: b.vel = spherical(velocity) * sqrt(scalefactor)
174: b
175:end
176:
177:
178:mkplummer(parse_command_line(options_text))
学生C：あー、コマンドラインオプションとかベクトル型とか、同じような感じですね。
赤木 ：もちろん、 Ruby のほうが元だから。
学生C：これならなんとかなるかな。やってみます。
赤木 ：お願いね。
学生C：一応できて動く気が、、、動かすと
gravity> crystal mkplummer.cr -- -n 5 -Y --- !CommandLog command: /home/makino/.cache/crystal/crystal-run-mkplummer.tmp -- -n 5 -Y log: Plummer model created --- !Particle id: 0 t: 0.0 m: 0.2 r: x: 0.23759233968153254 y: -0.3444540700300935 z: -0.48796274751905533 v: x: -0.5313327088566941 y: -0.11499495829365827 z: -0.4787121762076616 --- !Particle id: 1 t: 0.0 m: 0.2 r: x: -0.9682097409577199 y: 0.8313914547629354 z: 1.2489064796551024 v: x: 0.0852478489055863 y: -0.07597769258331093 z: 0.3888763253806409 --- !Particle id: 2 t: 0.0 m: 0.2 r: x: -0.1578932235094312 y: 0.038799916927366514 z: -0.2146892341645089 v: x: -0.17000567816595114 y: -0.2175635464933193 z: 0.23864512418885106 --- !Particle id: 3 t: 0.0 m: 0.2 r: x: 0.26009563389137885 y: -0.41452175988638007 z: -0.08577560770403372 v: x: 0.0659118990938399 y: -0.002376439304386317 z: -0.8460511179628695 --- !Particle id: 4 t: 0.0 m: 0.2 r: x: 0.6284149908942399 y: -0.1112155417738284 z: -0.4604788902675043 v: x: 0.5501786390232193 y: 0.41091263667467476 z: 0.6972418446010391こんな感じです。オプション2つあって、 -n で粒子数、 -Y で YAML でだす、 という感じです。プログラムは
1:require "clop"
2:require "./vector3.cr"
3:require "./nacsio.cr"
4:include Math
5:
6:optionstr = <<-END
7:
8: Description: Plummer's Model Builder
9: Long description:
10: This program creates an N-body realization of Plummer's Model.
11:
12: Original Ruby code:
13: (c) 2004, Piet Hut and Jun Makino; see ACS at www.artcompsi.org
14: The algorithm used is described in Aarseth, S., Henon, M., & Wielen, R.,
15: Astron. Astroph. 37, 183 (1974).
16:
17: Crystal Version (c) 2020- Jun Makino
18:
19: Short name:-n
20: Long name: --n_particles
21: Value type: int
22: Default value: 10
23: Variable name: n
24: Description: Number of particles
25: Long description:
26: Number of particles in a realization of Plummer's Model.
27: Each particles is drawn at random from the Plummer distribution,
28: and therefore there are no correlations between the particles.
29: Standard Units are used in which G = M = 1 and E = -1/4, where
30:
31: G is the gravitational constant
32: M is the total mass of the N-body system
33: E is the total energy of the N-body system
34:
35:
36: Short name: -s
37: Long name: --seed
38: Value type: int
39: Default value: 0
40: Variable name: seed
41: Description: pseudorandom number seed given
42: Long description:
43: Seed for the pseudorandom number generator. If a seed is given with
44: value zero, a preudorandom number is chosen as the value of the seed.
45: The seed value used is echoed separately from the seed value given,
46: to allow the possibility to repeat the creation of an N-body realization.
47:
48: Short name: -Y
49: Long name: --yaml_io
50: Value type: bool
51: Variable name: use_nacsio
52: Description: use nacs io format (yaml based)
53: Long description: use nacs io format (yaml based)
54:END
55:clop_init(__LINE__, __FILE__, __DIR__, "optionstr")
56:options=CLOP.new(optionstr,ARGV)
57:
58:module Nacsio
59:class Particle
60: def ekin # kinetic energy
61: 0.5*@mass*(@vel*@vel)
62: end
63: def epot(body_array) # potential energy
64: p = 0
65: body_array.each{ |b|
66: unless b == self
67: r = b.pos - @pos
68: p += -@mass*b.mass/sqrt(r*r)
69: end
70: }
71: p
72: end
73: def write
74: printf("%22.15e", @mass)
75: @pos.to_a.each { |x| printf("%23.15e", x) }
76: @vel.to_a.each { |x| printf("%23.15e", x) }
77: print "\n"
78: end
79:end
80:end
81:
82:include Nacsio
83:
84:class NBody
85:
86: property :time, :body
87: def initialize
88: @time = 0.0
89: @body = Array(Particle).new
90: end
91: def ekin # kinetic energy
92: e = 0
93: @body.each{|b| e += b.ekin}
94: e
95: end
96: def epot # potential energy
97: e = 0
98: @body.each{|b| e += b.epot(@body)}
99: e/2 # pairwise potentials were counted twice
100: end
101:
102: def adjust_center_of_mass
103: m_com = @body.reduce(0.0){|s, b| s + b.mass}
104: pos_com = @body.reduce(Vector3.new){|vec, b| vec + b.pos*b.mass}
105: vel_com = @body.reduce(Vector3.new){|vec, b| vec + b.vel*b.mass}
106: pos_com /= m_com
107: vel_com /= m_com
108: @body.each do |b|
109: b.pos -= pos_com
110: b.vel -= vel_com
111: end
112: end
113:
114: def adjust_units
115: alpha = -epot / 0.5
116: beta = ekin / 0.25
117: @body.each do |b|
118: b.pos *= alpha
119: b.vel /= sqrt(beta)
120: end
121: end
122: def write
123: print @body.size, "\n"
124: printf("%22.15e\n", @time)
125: @body.each do |b| b.write end
126: end
127: def nacswrite
128: @body.each{|b| print b.to_nacs}
129: end
130:
131:end
132:
133:def spherical(r, ran)
134: theta = acos(ran.rand(2.0)-1.0)
135: phi = ran.rand(2*PI)
136: Vector3.new( r * sin( theta ) * cos( phi ),
137: r * sin( theta ) * sin( phi ),
138: r * cos( theta ))
139:end
140:
141:def plummer_sample(r)
142: b = Particle.new
143: scalefactor = 16.0 / (3.0 * PI)
144: radius = 1.0 / sqrt( r.rand ** (-2.0/3.0) - 1.0)
145: b.pos = spherical(radius,r) / scalefactor
146: x = 0.0
147: y = 0.1
148: while y > x*x*(1.0-x*x)**3.5
149: x = r.rand
150: y = r.rand(0.1)
151: end
152: velocity = x * sqrt(2.0) * ( 1.0 + radius*radius)**(-0.25)
153: b.vel = spherical(velocity,r) * sqrt(scalefactor)
154: b
155:end
156:
157:if options.seed == 0
158: r= Random.new
159:else
160: r= Random.new(options.seed)
161:end
162:nb = NBody.new
163:options.n.times do |i|
164: b = plummer_sample(r)
165: b.mass = 1.0/options.n
166: b.id = i
167: nb.body.push(b)
168:end
169:nb.adjust_center_of_mass if options.n > 0
170:nb.adjust_units if options.n > 1 && options.n < 100000
171:if options.use_nacsio
172: print CommandLog.new("Plummer model created").to_nacs
173: nb.nacswrite
174:else
175: nb.write
176:end
177:
です。
赤木 ：解説一応お願いしていい？
学生C：もとの Ruby のを動くようにしただけなんであんまりわかってないで すが、一応やります。
56行目までは clop とか nacsio.cr とか読み込んで、あとコマンドラインオ プションの定義ですね。あ、 -s で乱数のシードを与える、というのもありま した。58行目から80行目までは Particle クラスを拡張して、運動エネルギー とポテンシャルエネルギーを計算させてます。ポテンシャルエネルギーは 粒子の配列とで全部計算ですね。この辺は、粒子の位置・速度のスケーリング に使ってます。
で、粒子系のクラスというのを作ってあって、それが84行目からですね。 これは全体の運動エネルギー、ポテンシャルエネルギーを計算する 関数 ekin, epot とか、重心速度、重心位置を原点にもってくる adjust_pos、ポテンシャルエネルギーと運動エネルギーをそれぞれ -0.5, 0.25 になるようにスケールする ajust_units と、1行1粒子とかYAMLとかで書 くwrite, nacswrite ですう。中身はまあみての通りで。
spherical は球面上の一様乱数ですかね。ここで、 rand は Crystal の乱数 で、引数に実数を与えると0からそこまでの乱数になります。
plummer_sample が実際にプラマーモデルの分配に従って粒子1つを生成する関 数です。あんまりわかってないですが論文通りなんじゃないかと、、、元の Ruby のプログラムからあんまりここいじってなくて、乱数の関数の書き方か えたくらいなので。
あとは、157行からの if で乱数の初期化をして、163行からで n 個粒子作っ て、あと重心を0にとかエネルギーのスケーリングとかして出力、という感じ ですね。
赤木 ：これ、粒子数増やしてプロットできる？
学生C：オプションでですか？
赤木 ：あ、じゃなくて、折角だからこの出力ファイルを読んでプロットするのにして。プロットするの自体は前にやってるわよね」？
学生C：ですね。読むのもさっきやったから、、、こんなのですね。
Figure 26: 1000粒子プラマーモデルの出力
crystal mkplummer.cr -- -n 1000 -Y | nacsplotで作ってて、nacsplot は
1:require "grlib"
2:require "clop"
3:require "./nacsio.cr"
4:include Math
5:include GR
6:include Nacsio
7:
8:optionstr= <<-END
9: Description: Plot program for nacs snapshot
10: Long description: Plot program for nacs snapshot
11:
12: Short name:-w
13: Long name: --window-size
14: Value type: float
15: Variable name:wsize
16: Default value:1.5
17: Description:Window size for plotting
18: Long description:
19: Window size for plotting orbit. Window is [-wsize, wsize] for both of
20: x and y coordinates
21:END
22:
23:clop_init(__LINE__, __FILE__, __DIR__, "optionstr")
24:options=CLOP.new(optionstr,ARGV)
25:update_commandlog
26:pp=Array(Particle).new
27:while (sp= CP(Particle).read_particle).y != nil
28: pp.push sp.p
29:end
30:wsize=options.wsize
31:setwindow(-wsize, wsize,-wsize, wsize)
32:setcharheight(0.05)
33:box
34:mathtex(0.5, 0.06, "x")
35:mathtex(0.06, 0.5, "y")
36:text(0.6,0.91,"t="+sprintf("%.3f",pp[0].time))
37:setmarkertype(4)
38:setmarkersize(1)
39:polymarker(pp.map{|p| p.pos[0]}, pp.map{|p| p.pos[1]})
40:c=STDERR.gets
です。説明いります？だいたいみての通りですが、、、
赤木 ：お願いね。
学生C：24行目までは色々なライブラリとかと、あとコマンドラインオプショ ンです。ここでは -w でプロットする座標範囲、というのだけ残してます。
で、CommandLog 読まないといけないので update_commandlog を読んで、 それから26行目から29行目で粒子の配列を作ります。このプログラムは粒子デー タ出力しないので、Particle クラスのほうだけとって配列にいれます。
30行目からは前の particle1.cr かなんかと同じです。あ、 最後の、入力待つところ、標準入力は粒子データのリダイレクトにしているの で、それでもキーボードの入力待ちになるように STDERR にしてます。
赤木 ：こういうのがあると、 x, y 座標に色々指定できるようにしたくなる わね。
学生C：コマンドラインで式を与えるとかですか？そういうのは Ruby とか Python でしたほうがいいものではないですか？
赤木 ：でもまあ100倍の速度の違いは問題で、そうするとやっぱり大きなデー タ使うと遅いとかになっちゃうから、、、
学生C：ちょっと話が長くなりそうなので次回以降にしませんか？
赤木 ：そうね。
13.1. 課題
累積密度分布は、半径方向に粒子をソートして、各粒子の位置までで内側にあ
る質量を求めればよい。ソートには Array#sort! を使う。 sortを使うには、
演算子 <=> を定義するか、あるいはそれに相当するブロックを渡す。以下は
ブロックを渡す例である。
data= [9,6,7,5,4,8,2,3,0,1]
data.sort!{|x,y|
if x>y
1
elsif x<y
-1
else
0
end }
p! data
等しいなら0、大小に応じて 負または正の値を返すようにする。
13.2. まとめ
13.3. 参考資料
| Previous | ToC | Next |