概要
初投稿です。卒業研究でMATLABを用いたモンテカルロシミュレーションを行っていました。そこで任意の形状内を散乱する電子の経路を求めるプログラムを作成したので、まとめたいと思います。
モンテカルロ法
モンテカルロ法とは、簡単に、乱数を活用したシミュレーションの手法です。今回の場合では、電子の散乱方向やタイミングに乱数を用いています。電子がたくさんある場所では、電子同士がぶつかり散乱をします。電子の経路を求めるためには他の大量の電子の動きが影響するため、計算量が莫大になります。そこで乱数を活用することで、電子の散乱の経路を他の電子の影響を受けずに計算することができます。
実装
ここから実際の実装したプログラムを用いて説明していきます。
電子数は10~100万個程で計算していました。経路の計算に他の電子は影響しないので、parforを用いて、GPUを活用した並列計算を行っています。それでは説明していきます。
形状設計
形状は直線の頂点を、記録しています。polyshape関数で形を作った後、boundary関数でその頂点を配列に記録します。(ape_x, ape_y) また後に必要なため、傾きも記録しています。(slp)
%% 正多角形 (円、楕円、正方形 etc.)
th = 60:60:360;
pgon = polyshape(x10e-6*cos(deg2rad(th)),10e-6*sin(deg2rad(th)));
[ape_x, ape_y] = boundary(pgon); % 壁の頂点の座標を配列に記録
slp = zeros(1, length(ape_x) - 1); % 壁の傾きを配列に記録
for i = 1:length(ape_x) - 1
slp(i) = atan((ape_y(i + 1) - ape_y(i)) /(ape_x(i + 1) - ape_x(i)));
end
構造内に電子を用意
自分の研究では、光を当てた電子の経路を計算します。そのため光の当たった電子を構造内に入れます。光の当たる電子は、光の中心から正規分布になります。現在入っている電子の数を記録しておき、目標の個数入るまで、正規分布した乱数で電子の座標を決めます。
- inpt_pos_x,inpt_pos_y : 光を当てる中心の座標が入った配列。
- inpt_num : 構造内にある電子の個数です。
- pump : 光のサイズで正規分布の分散です。
- x_cand,y_cand : 正規分布した乱数で決まった電子の座標です。
- one_pos_ele : 1か所当たりの電子の目標個数です。
for i = 1 : length(inpt_pos_x) % 光の場所それぞれ行います。
while inpl_num < i*one_pos_ele %% 目標個数入るまで続けます。
x_cand = inpt_pos_x(i) + pump.*randn(1,i*one_pos_ele - inpl_num);
y_cand = inpt_pos_y(i) + pump.*randn(1,i*one_pos_ele - inpl_num);
[in,~] = inpolygon(x_cand,y_cand,ape_x,ape_y);
% 構造内だった電子だけを抜き取ります。
for j = 1:sum(in) % 構造内だったものを追加します。
x(inpl_num+1 : inpl_num+sum(in)) = x_cand(in~=0);
y(inpl_num+1 : inpl_num+sum(in)) = y_cand(in~=0);
end
inpl_num = inpl_num + sum(in); % 構造内の電子数を更新します。
end
end
散乱経路を求める
ある一つの電子の経路を求めます。
最終的に知りたいのは、特定の時間における電子のx座標とy座標です。電子は、散乱と壁での反射によって方向を変えながら基本的に等速直進運動をすると考えます。そのため、電子が方向を変える座標、時間を計算し、その後、特定の時間の座標を求めます。
- To_next_scattering() : 方向を変える座標、時間を記録していく関数です。
- Chg_ang_X, .._Y, .._time : 方向が変わった時の座標、時間が入った配列です。
- Cnt_chg_ang : 方向が変わった回数、上の3つの配列の大きさです。
- scattering_angle : 散乱後の角度が入った配列、2π x 0~1 の乱数
- scattering_interval : 散乱-散乱間の時間が入った配列
- nowtime : 現在時刻
引数が-1になっているのは、反射したときに壁の番号を記録するものです。散乱後や最初は-1にしています。
for j = 1:length(scattering_point)
[X(j+1), Y(j+1), Chg_ang_X, Chg_ang_Y, Chg_ang_time, nowtime, Cnt_chg_ang] = ...
To_next_scattering(X(j), Y(j), scattering_angle(j), scattering_interval(j), vF, ape_x, ape_y, slp, ...
Chg_ang_X, Chg_ang_Y, Chg_ang_time, nowtime, Cnt_chg_ang, -1);
end
To_next_scattering関数は次のようになっています。
まず、仮に壁がなかったときの次の散乱の座標(x_bal,y_bal)を、求め壁一つ一つとそこに到達するまでにぶつかるかをCrossing関数で判別します。そのbool値をany_reflectという配列に入れます。
ps = 1e-12;
x_bal = x + cos(angle)*interval*vF*ps;
y_bal = y + sin(angle)*interval*vF*ps;
any_reflect = true(1, length(ape_x) - 1);
for i = 1:length(any_reflect)
if(i == num)
any_reflect(i) = false;
else
[any_reflect(i)] = Crossing(x, y, x_bal, y_bal, ape_x(i), ape_y(i), ape_x(i+1), ape_y(i+1), any_reflect(i));
end
end
- angle : 電子の進む方向
- interval : 次の散乱までの時間
Crossing関数は次のようになります。線分abと線分cdの交差判定は、簡単な数学で解けます。s,tはそれぞれ、cとdがabに対してどっち側にいるかを表します。sとtの符号が異なることは、別側にいる、つまり交差してることを表します。同じことをa,bとcdで行い、どちらかが別側なら交差しているとしてjudge=trueを返します。
function [judge] = Crossing(ax,ay,bx,by,cx,cy,dx,dy,judge)
s = (ax - bx) * (cy - ay) - (ay - by) * (cx - ax);
t = (ax - bx) * (dy - ay) - (ay - by) * (dx - ax);
if (s * t > 0)
judge = false;
end
s2 = (cx - dx) * (ay - cy) - (cy - dy) * (ax - cx);
t2 = (cx - dx) * (by - cy) - (cy - dy) * (bx - cx);
if (s2 * t2 > 0)
judge = false;
end
end
交差し得る壁がわかりました。それが存在するとき、最初にぶつかる壁で反射します。その壁を見つけていきます。
まず交差し得る壁との交点を求めます(intersec_x,intersec_y)。そしてその距離を求め(dist)、それが一番近い壁ならば、その距離、壁の番号、壁のどこで反射するかを記録します(far_from_first_wall, nearest_wall, ref_x, ref_y)。
nearest_wall = -1;
far_from_first_wall = 1000000;
for i = 1:length(any_reflect)
if(any_reflect(i))
[intersec_x, intersec_y] = Get_intersect(x,y,x_bal,y_bal,ape_x(i),ape_y(i),ape_x(i+1),ape_y(i+1));
dist = DistanceBet2(x,y,intersec_x,intersec_y);
if(far_from_first_wall > dist)
far_from_first_wall = dist;
nearest_wall = i;
ref_x = intersec_x;
ref_y = intersec_y;
end
end
end
ここで交点を求めるGet_intersect関数、距離を求めるDistanceBet2関数はただの簡単な数学です。
function [intersec_x, intersec_y] = Get_intersect(x,y,x2,y2,x3,y3,x4,y4)
det = (x - x2) * (y4 - y3) - (x4 - x3) * (y - y2);
t = ((y4 - y3) * (x4 - x2) + (x3 - x4) * (y4 - y2)) / det;
intersec_x = t * x + (1.0 - t) * x2;
intersec_y = t * y + (1.0 - t) * y2;
end
function [distance] = DistanceBet2(x,y,x2,y2)
distance = sqrt((x2-x)*(x2-x) + (y2-y)*(y2-y));
end
ここまでで次の散乱までに反射が起こるか、起こるなら反射する壁はどれか、がわかりました。起こらないなら、関数を終わります。そして、起こるなら、反射する場所を始点としてもう一度、To_next_scattering関数を呼びます。反射した壁の番号(nearest_wall)を引数に入れ、次はその壁で反射しないようにします。
cnt = cnt + 1; % 方向を変更した合計回数を更新します。
if(nearest_wall == -1) % 反射が起こらない場合です。
% 次に電子の方向が変わる座標、時間を更新します
x_next = x_bal;
y_next = y_bal;
nowtime = nowtime + interval;
% 方向が変わった座標、時間を配列に追加します。
chg_x(cnt) = x_bal;
chg_y(cnt) = y_bal;
chg_time(cnt) = nowtime;
else % 反射が起こる場合です。
% 反射後の電子の進む方向を求めます。
ref_angle = mod(2*slp(nearest_wall) - angle, 2*pi);
% 反射する座標までにかかる時間を求めます。
elp_time = (far_from_first_wall / DistanceBet2(x, y, x_bal, y_bal)) * interval;
% 反射する時間で更新します。
nowtime = nowtime + elp_time;
% 方向が変わる座標、時間を配列に追加します。
chg_x(cnt) = ref_x;
chg_y(cnt) = ref_y;
chg_time(cnt) = nowtime;
% 反射点を始点として、もう一度関数を呼びます。
[x_next, y_next, chg_x, chg_y, chg_time, nowtime, cnt] = ...
To_next_scattering(ref_x, ref_y, ref_angle, interval - elp_time, vF, ape_x, ape_y, slp, chg_x, chg_y, chg_time, nowtime, cnt, nearest_wall);
end
観測終了時間まで、To_next_scattering関数で、Chg_ang_X,Chg_ang_Y,Chg_ang_timeを更新していきます。それによって、電子の方向が変わったときの、座標と時間は分かり。そして、その間は直線運動をするため、電子の経路がわかり、任意の時間の電子の座標がわかります。下の画像は、反射が正しくできてるかを表すものです。円(36万角形)内を反射しながら進み、中心にまた円ができています。
これに散乱を加えると下のようになります。正六角形内で散乱と反射を繰り返します。
最後に
このようにして、任意の形状内を散乱する電子の経路が求められました。
できるだけ、簡単な数学的な問題にして、求められるように意識しました。