3次元井戸型ポテンシャルのパイロットウェーブ、つくりかた

これの作り方備忘録です。


図のように、ｘ、ｙ、ｚ、ｒの列を作り、そこに一様乱数を叩き込みます。
今回は1000点ほどデータを取りました。
ｒの列はrand()なので、0～1の一様乱数ですが
ｘ、ｙ、ｚは-1～1の一様乱数にしたかったので、2*rand()-1を使ってます。


次に、井戸型ポテンシャルにおける波動関数、パイロットウェーブのようなものを作ります。

左3列がそれぞれｘ、ｙ、ｚによる偶関数の波動関数
同様に右3列がｘ、ｙ、ｚによる奇関数の波動関数です。

偶関数の場合は
cos(2πnx/4)を
奇関数の場合は
sin(2πnx/4)を、
それぞれ使っています。ｙ、ｚの波動関数(H,I,K,L列)は、xの代わりにyやzの乱数座標を入れてください。

また、ここでいうnとは、4行目にある数字のことを言いますが
1行目に入れた数が偶数ならG,H,I列をゼロに、奇数ならJ,K,L列をゼロにするように
2を法としたモジュロ演算を使っています。

波動関数そのものも、
1行目にいれたx,y,zの量子数が偶数ならG,H,I列の波動関数をゼロに、奇数ならJ,K,Lをゼロにするようにしています。

ここでは、x軸とz軸の量子数が2と4と、偶数なので
G,I列がゼロになって
K,L列に有限の数を入れました。

逆に、y軸の量子数が3と奇数なので
H列に有限の値、I列にゼロが入っています。

 

次に、各軸の波動関数を完成させます。
 
偶関数か奇関数のどちらかが入るはずなので、
M列=G列+J列
としています。
y,z軸もそれぞれ、N列=H列+K列、O列=I列+L列としています。


それから、3次元の波動関数を完成させたいので
一旦変数分離した各軸の波動関数を掛け合わせます。
P列=M列*N列*O列といった感じです。Ψ=Ψx*Ψy*Ψzですね。

確率は、Ψの2乗にrのパイロットウェーブ的なのを掛け算したいので
Q列=P列^2*F列としたいところですが、
今回の場合規格化が無理そうなので、代わりに、波動関数に定数を掛け算して、ちょうどよくなるように調整しました。
Q4に適切な値を入れて、掛け算しているので

Q列=P列^2*F列*$Q$4
となります。

最後に、R列で、Q列の値が0.5未満だったら0を、0.5以上だったら1を返すようにしました。
敷居値をR4に設けつつ、if文を使用しましたが、
別に0桁でのround関数で構いません。

この1の個数が、値の個数の半分くらいになるように、Q4の値を調整しています。
結果的に40ぐらいが適切と判断しました。




このシートのすべてをコピーし、新しいシートに、数式ではなく値として貼り付けます。


ｘ、ｙ、ｚ、ｐの値だけを使うのですが、ｐが1のときのデータだけがほしいので、ｐの列に関してデータを降順に並べ替えます。
並べ替えるには、この4列をひとまとまりのデータと機械に認識してもらうために
周りを空けて、ｘ、ｙ、ｚ、ｐの列ラベルに2つ以上の書式変更を設定します。
今回は中央寄せと太文字を用いました。
また、降順などの操作を行う際に、数式コピペだとうまくいかないことが多いので、一旦値コピペをしました。(ほかにも理由はあります)

これで、ｐが1の区間の行すべてを取り囲み、グラフ化することで、外村彰先生の実験のような
パイロットウェーブ的な干渉縞をシミュレーションすることができます。

ただし、ここで用いるグラフは2Dグラフで、ｘ、ｙのデータしか用いません。
ですので、ｘ、ｙ、ｚのデータを適切に回転させることで、3Dのように見せかけることにしました。

図のような、ｙ軸回転のあとにｘ軸回転させるような回転行列を作用させます。

H8～J10の区間に、mmult(H2:J4,H5:J7)という行列の積を計算させます。
まず代表となるH8セルに行列の積を実装し、H8～J10を選択して
ｃｔｒｌ＋ｓｈｉｆｔ＋Ｅｎｔｅｒを押すことで、H8～J10全体に計算結果をいきわたらせることができます。

そしてこのH8～J10の混合された回転行列を、先ほどのC,D,E列に作用させたいので
mmult(C12:E12,$H$8:$J$10)とします。
C～E列の500列くらいに作用させるのに、H8～J10の行列は固定のまま作用させたいので
ここでも絶対参照を用います。


これで、G3やG6セルをいじれば、ｘ、ｙ平面だけでなく、見たい角度からの粒子を見ることができるようになりました。

このG3やG6セルも、自動で動かしたいため、以下のような表を作ります。

P列の値に応じて、Q列のｙ軸回転の角度(°)と、R列のｘ軸回転角度(°)をまとめてあります。
Q3とR3セルには、それぞれ、選択したP列の数字に応じたQ列とR列の値が出るように
vlookup関数を用いています。

ここで、選択したP列の数字を時々刻々と動かしたいので
now()-today()を計算して、シリアル値を0～1までの数値として表現します。
この値は隠れがちですが1秒以内に何度も更新されています。
これを100万倍ぐらいさせて、わずかな値の変化を増幅します。
それからround関数で整数にしてやり
先ほどの図のP列の最大値を法としたモジュロ演算を行うことで、周期をP列の値の整数にすることができます。ここでは89行あったので89周期です。

ただし、mod関数は0から88までの89個が算出されるので、1を足して、1から89までの整数に直してやります。

この結果が仮にC6セルにあったとすると
Q3とR3にはそれぞれ
vlookup(C6,P5:R93,2)
vlookup(C6,P5:R93,3)
と入力することで、C6セルがP列の値と等しい行のQ列とR列の値を算出してくれます。
(vlookupの3つ目の変数である2や3は、参照する範囲の2列目、3列目という意味です)

あとは
HとI列を参照して、散布図のグラフにすれば完成です


Excelファイルをアップしておきます。
windowsのPCでしか動作確認してませんが、テキトーな空白セルで、delボタンを連打していただければそのたびに再計算されるので、グラフが動くと思います