tobiom’s diary

備忘録

gnuplotを使ったCoulomb散乱の図示化

gnuplotでCoulomb散乱のモデルを計算してコマ毎に描画したのをWindowsムービーメーカーで繋げて動画にしてみました。
計算には2次のSymplectic法を用いています。

コードを書いて

Coulomb_Scattering.plt
#Coulomb_Scattering

set terminal jpeg enhanced font "Times" 20 size 600, 480
set tics font 'Times,12'
set nokey
set ticslevel 0

#画面
wx=10.0
wy=8.0
set xr[-wx*0.5:wx*0.5]
set yr[-wy*0.5:wy*0.5]

ASIZE=2*15

#配列にデータを格納する関数を定義するための文字列の生成
str_qx = "qx_set(i,x) = "
str_qy = "qy_set(i,x) = "
str_px = "px_set(i,x) = "
str_py = "py_set(i,x) = "
do for [k=0:ASIZE-1] {str_qx = str_qx . sprintf("(int(i)==%d ? _qx%d=x :", k, k)}
do for [k=0:ASIZE-1] {str_qy = str_qy . sprintf("(int(i)==%d ? _qy%d=x :", k, k)}
do for [k=0:ASIZE-1] {str_px = str_px . sprintf("(int(i)==%d ? _px%d=x :", k, k)}
do for [k=0:ASIZE-1] {str_py = str_py . sprintf("(int(i)==%d ? _py%d=x :", k, k)}
str_qx = str_qx. "1/0"
str_qy = str_qy. "1/0"
str_px = str_px. "1/0"
str_py = str_py. "1/0"
do for [k=0:ASIZE-1] {str_qx = str_qx . ")"}
do for [k=0:ASIZE-1] {str_qy = str_qy . ")"}
do for [k=0:ASIZE-1] {str_px = str_px . ")"}
do for [k=0:ASIZE-1] {str_py = str_py . ")"}
#文字列実行
eval str_qx
eval str_qy
eval str_px
eval str_py

#配列のデータを参照する関数を定義するための文字列の生成
str_qx = "qx_get(i) = "
str_qy = "qy_get(i) = "
str_px = "px_get(i) = "
str_py = "py_get(i) = "
do for [k=0:ASIZE-1] {str_qx = str_qx . sprintf("(int(i)==%d ? _qx%d :", k, k)}
do for [k=0:ASIZE-1] {str_qy = str_qy . sprintf("(int(i)==%d ? _qy%d :", k, k)}
do for [k=0:ASIZE-1] {str_px = str_px . sprintf("(int(i)==%d ? _px%d :", k, k)}
do for [k=0:ASIZE-1] {str_py = str_py . sprintf("(int(i)==%d ? _py%d :", k, k)}
str_qx = str_qx. "1/0"
str_qy = str_qy. "1/0"
str_px = str_px. "1/0"
str_py = str_py. "1/0"
do for [k=0:ASIZE-1] {str_qx = str_qx . ")"}
do for [k=0:ASIZE-1] {str_qy = str_qy . ")"}
do for [k=0:ASIZE-1] {str_px = str_px . ")"}
do for [k=0:ASIZE-1] {str_py = str_py . ")"}
#文字列実行
eval str_qx
eval str_qy
eval str_px
eval str_py

DATA = "Coulomb_Scattering.data"
set print DATA

#初期値・係数
mass = 0.5
alpha = 1.5
R=wx*0.4
v0 = 3.0
do for [k=0:ASIZE-1] {dummy=qx_set(k,-wx*0.5)}
do for [k=0:ASIZE-1] {dummy=qy_set(k,(k-(ASIZE-1)*0.5)*0.2)}
do for [k=0:ASIZE-1] {dummy=px_set(k,mass*v0)}
do for [k=0:ASIZE-1] {dummy=py_set(k,0)}

#微分方程式
partial_p(p) = (p)/(mass)
partial_qx(qx,qy) = -alpha*(qx-R)/(((qx-R)**2+qy**2)**(3/2))
partial_qy(qx,qy) = -alpha*qy/(((qx-R)**2+qy**2)**(3/2))
dt = 0.001    #時間刻み
n_max = 500  #静止画数
m_max = 20   #間引き数
n = 0;
m = 0
call "si2.plt"
si2.plt
do for [k=0:ASIZE-1] {
	cpx=px_get(k)
	cpy=py_get(k)
	cqx=qx_get(k)
	cqy=qy_get(k)
	dummy=px_set(k,cpx-dt*0.5*partial_qx(cqx,cqy))
	dummy=py_set(k,cpy-dt*0.5*partial_qy(cqx,cqy))
}
do for [k=0:ASIZE-1] {
	cpx=px_get(k)
	cpy=py_get(k)
	cqx=qx_get(k)
	cqy=qy_get(k)
	dummy=qx_set(k,cqx+dt*partial_p(cpx))
	dummy=qy_set(k,cqy+dt*partial_p(cpy))
}
do for [k=0:ASIZE-1] {
	cpx=px_get(k)
	cpy=py_get(k)
	cqx=qx_get(k)
	cqy=qy_get(k)
	dummy=px_set(k,cpx-dt*0.5*partial_qx(cqx,cqy))
	dummy=py_set(k,cpy-dt*0.5*partial_qy(cqx,cqy))
}

do for [k=0:ASIZE-1] {print qx_get(k),qy_get(k)}
m=m+1
if(m<m_max) reread
m=0

outfile(n) = sprintf("f%d.jpg",n+10000) #出力ファイル名
title(n) = sprintf("t = %4.2f",n*dt) #タイトル名
unset label
set label title(n) font 'Times,12' at -wx/2, wy/2
set output outfile(n)
plot DATA with points pt 1 ps 0.1

n = n+1
if(n < n_max) reread
n=0

実行して画像を生成

f:id:tobiom:20161216234524p:plain
f:id:tobiom:20161216234445j:plain

参考文献

gnuplotでの配列(もどき)の使用 - 米澤進吾 ホームページ
gnuplotには配列が標準装備で無かったのでこのサイトのコードをコピペして改悪参考にさせて頂きました。