球面調和関数（４）

 球面調和関数（３）からの続き

多くの疑問はよく探してみると答えになるウェブサイトがみつかるものだ。Juliaによる球面調和関数というか原子波動関数の可視化のための関数定義もここにある。面倒なのがラゲール陪関数の定義や規格化なのだけれど，とりあえずなんとかなった。あとは球面調和関数で磁気量子数が負の場合の扱いなど。

using GSL
using Plots
using LinearAlgebra

function Y(l,m,θ,ϕ)
# spherical harmonics
# l: orbital qn, m: magnetic qn
# θ,φ: spherical angular coordinate
#
  am=abs(m)
  if m <= 0 || mod(am,2) == 0
    p=1
  else
    p=-1
  end
#
  if am > l
    return 0
  else
    return p*sf_legendre_sphPlm(l,am,cos(θ))*exp(im*am*ϕ)
  end
end

function R(n,l,r)
# Z: nuclear charge (set 1 for hydrogen)
# n: principle qn, l: orbital qn
# r: Bohr radius unit a0 = ℏc/α・mc^2
#
Z=1
ρ=2*Z*r/n
  if n < l
    return 0
  else
    return sf_laguerre_n(n-l-1,2*l+1,ρ)*exp(-ρ/2)*ρ^l
  end
end

function norm(n,l)
# square of radial wave function norm
  return (2/n)^3 *factorial(n-l-1)/(2*n*factorial(n+l))
end

x=(0.0:0.1:25)
#
r10=(R.(1,0,x).*x).^2 .*norm(1,0)
r20=(R.(2,0,x).*x).^2 .*norm(2,0)
r21=(R.(2,1,x).*x).^2 .*norm(2,1)
r30=(R.(3,0,x).*x).^2 .*norm(3,0)
r31=(R.(3,1,x).*x).^2 .*norm(3,1)
r32=(R.(3,2,x).*x).^2 .*norm(3,2)
#
plot(x,r10,linewidth=1,legend=:topright)
plot!(x,[r20,r21],linewidth=1.5)
g1=plot!(x,[r30,r31,r32],linewidth=2,linestyle=:dash)

t1=(0.0:0.01:pi)
t2=(pi:0.01:2pi)
#
y10=abs2.(Y.(1,0,t1,0))
y11=abs2.(Y.(1,1,t1,0))
y20=abs2.(Y.(2,0,t2,0))
y21=abs2.(Y.(2,1,t2,0))
y22=abs2.(Y.(2,2,t2,0))
#
plot(t1,[y10,y11],proj=:polar,linewidth=2)
g2=plot!(t2,[y20,y21,y22],proj=:polar,linewidth=2)

plot(g1,g2)

図　水素原子の動径部分と角度部分