球面調和関数の可視化【SageMath】【Python】

球面調和関数 $Y_l^m(\phi,\theta) \ \ (0\leq \phi <2\pi,0\leq \theta <\pi)$ の可視化を行います。量子力学でいうと l は方位量子数、m は磁気量子数になります。

用いる式は以下のとおりです。

$r(\phi,\theta)=|Y_l^m(\phi,\theta)|^2\\ x=r\sin(\theta)\cos(\phi)\\ y=r\sin(\theta)\sin(\phi)\\ z=r\cos(\theta)$

自明ですが、r が定数の場合は球面となります。

l=3,m=1

SageShell を起動し、以下のとおり入力します。

sage: phi,theta = var('phi,theta')
sage: Y31 = spherical_harmonic(3,1,theta,phi)
sage: r = (Y31(phi,theta).conjugate()*Y31(phi,theta)).real()
sage: parametric_plot3d((r(phi,theta)*sin(theta)*cos(phi),r(phi,theta)*sin(theta)*sin(phi),r(phi,theta)*cos(theta)),(phi,0,2*pi),(theta,0,pi))
Launched html viewer for Graphics3d Object

結果はブラウザで表示されます。

f:id:riverta1992:20210102181954p:plain

$Y_l^m$ は複素数値をとる関数なので $r(\phi,\theta)=|Y_l^m(\phi,\theta)|^2$ では複素数の重要な情報が抜け落ちてしまいます。ここでは複素数の偏角と曲面の色を対応させることとし、再度可視化を行います。

まず、以下で複素数 n の偏角を求める関数を定義します。

sage: def my_arg(n):
....:     rad = 0.0
....:     if bool(n.real()==0):
....:         rad = pi/2
....:     elif bool(n.imag()>0) and bool(n.real()>0):
....:         rad = arctan(n.imag()/n.real())
....:     elif bool(n.imag()>0) and bool(n.real()<0):
....:         rad = arctan(n.imag()/n.real())+pi
....:     elif bool(n.imag()<0) and bool(n.real()>0):
....:         rad = arctan(n.imag()/n.real())+2*pi
....:     elif bool(n.imag()<0) and bool(n.real()<0):
....:         rad = arctan(n.imag()/n.real())+pi
....:     return float(rad/2/pi)
....:

偏角と色を対応づけるためにHSV色空間を用いて、色を指定する関数を定義します。

sage: cm = colormaps.hsv
sage: def c(phi,theta):return my_arg(Y31(phi,theta))

結果をプロットします。

sage: parametric_plot3d((r(phi,theta)*sin(theta)*cos(phi),r(phi,theta)*sin(theta)*sin(phi),r(phi,theta)*cos(theta)),(phi,0,2*pi),(theta,0,pi),color=(c,cm))

f:id:riverta1992:20210102174217p:plain

l=3,m=0

sage: Y30 = spherical_harmonic(3,0,theta,phi)
sage: r = (Y30(phi,theta).conjugate()*Y30(phi,theta)).real()
ValueError: the number of arguments must be less than or equal to 1

r の定義でエラーが出ます。Y30 の引数の数で怒られているようです。

sage: Y30.arguments()
(theta,)

Y30 がシータだけの関数になっていました。単にY30と入力するだけでも確認できます。

sage: Y30
-1/4*(5*sqrt(7)*cos(theta)*sin(theta)^2 - 2*sqrt(7)*cos(theta))/sqrt(pi)

これに注意してもう一度 r を定義します。

sage: r = (Y30(theta).conjugate()*Y30(theta)).real()

色と偏角の対応づけも再度行います。

sage: def c(phi,theta):return my_arg(Y30(theta))

最後にプロットします。

sage: parametric_plot3d((r(theta)*sin(theta)*cos(phi),r(theta)*sin(theta)*sin(phi),r(theta)*cos(theta)),(phi,0,2*pi),(t
....: heta,0,pi),color=(c,cm))
/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/repl/rich_output/display_manager.py:586: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...)
See http://trac.sagemath.org/5930 for details.
  return obj._rich_repr_(self)
Launched html viewer for Graphics3d Object

f:id:riverta1992:20210102192039p:plain

色が曲面のどの点も同じになっています。Y30 の値域が実数になっているので、どの点でも偏角が同じになるからですね。

参考

仮想計算機構

IT業界と無縁な派遣社員のブログ

球面調和関数の可視化【SageMath】【Python】

l=3,m=1

l=3,m=0

参考

球面調和関数

Sage関係

その他