原子f軌道のリスト

3次元グラフの表面の色づけの項で説明した「光の当たっている物体に見える色のつけ方」を利用して、原子軌道のうちf軌道の角度部分Y(theta,phi)のリストを作成しました。下のスクリプトを実行すると、こちら（4MBあります）のような4ページのPostScriptファイルが出来ます。それぞれのページに上図のように色使いの違う図が出来ます。

なお、このスクリプトの実行にはgnuplot 4.3以降が必要です。また、以下のスクリプトは図の解像度を上げているため実行にはかなり時間がかかります。解像度を下げて実行速度を上げたい場合はset isosamples 90とset samples 45の部分の数字を小さくしてください。 

#version 4.3WN09Sep

set encoding iso_8859_1

if (exist("n")==0) n=0

if (n==0 || n==1) set palette defined ( 0 "dark-blue" , 1 "cyan")
if (n==2) set palette defined ( 0 "brown" , 1 "yellow")
if (n==3) set palette defined ( 0 "dark-magenta" , 1 "pink")
if (n==4) set palette defined ( 0 "dark-green" , 1 "light-green")

h=0.000001

lx=2.0
ly=1.0
lz=1.0

# 表面の媒介変数関数
Fx(u,v)=sin(u)*cos(v)*abs(Y(u,v))
Fy(u,v)=sin(u)*sin(v)*abs(Y(u,v))
Fz(u,v)=cos(u)*abs(Y(u,v))

# 光の方向
sx=-1.0
sy=1.0
sz=1.0

# 微分用の微小数
h=0.000001

# Fx,Fy,Fzのuおよびvによる偏微分
dFxu(u,v)=(Fx(u+h,v)-Fx(u-h,v))/(2.0*h)
dFyu(u,v)=(Fy(u+h,v)-Fy(u-h,v))/(2.0*h)
dFzu(u,v)=(Fz(u+h,v)-Fz(u-h,v))/(2.0*h)
dFxv(u,v)=(Fx(u,v+h)-Fx(u,v-h))/(2.0*h)
dFyv(u,v)=(Fy(u,v+h)-Fy(u,v-h))/(2.0*h)
dFzv(u,v)=(Fz(u,v+h)-Fz(u,v-h))/(2.0*h)

# 法線ベクトルの各成分 n = lu x lv
nx(u,v)=dFyu(u,v)*dFzv(u,v)-dFyv(u,v)*dFzu(u,v)
ny(u,v)=dFzu(u,v)*dFxv(u,v)-dFzv(u,v)*dFxu(u,v)
nz(u,v)=dFxu(u,v)*dFyv(u,v)-dFxv(u,v)*dFyu(u,v)

# 光の方向ベクトルsの長さ
s_length=sqrt(sx**2+sy**2+sz**2)
# 法線ベクトルnの長さ
n_length(u,v)=sqrt(nx(u,v)**2+ny(u,v)**2+nz(u,v)**2)

# 光の方向ベクトルsと表面の法線ベクトルnの内積で色を決める
light(u,v)=(n_length(u,v)!=0.0 ? (sx*nx(u,v)+sy*ny(u,v)+sz*nz(u,v))/(s_length*n_length(u,v)) : light(u+h,v+h))

# 球面調和関数の実部と虚部
Y30(u,v)=0.25*sqrt(7.0/pi)*(5.0*cos(u)**3-3.0*cos(u))
Y31_re(u,v)=(-0.125)*sqrt(2.0*21.0/pi)*sin(u)*(5.0*cos(u)**2-1.0)*cos(v)
Y31_im(u,v)=(-0.125)*sqrt(2.0*21.0/pi)*sin(u)*(5.0*cos(u)**2-1.0)*sin(v)
Y32_re(u,v)=0.25*sqrt(105.0/pi)*sin(u)**2*cos(u)*cos(2.0*v)
Y32_im(u,v)=0.25*sqrt(105.0/pi)*sin(u)**2*cos(u)*sin(2.0*v)
Y33_re(u,v)=(-0.125)*sqrt(2.0*35.0/pi)*sin(u)**3*cos(3.0*v)
Y33_im(u,v)=(-0.125)*sqrt(2.0*35.0/pi)*sin(u)**3*sin(3.0*v)


set parametric
set isosamples 90
set samples 45
set angles radians

# 一列目にtheta、二列目にphiの値の入ったファイルを作る。
if (n==0) tablefile="theta-phi.dat";\
     set urange [0:pi];\
     set vrange [0:2*pi];\
     set table tablefile;\
     splot u,v,0;\
     unset table


set multiplot #layout 4,2 rowsfirst scale 1.6,1.6 offset 0.05,0.08


set object 1 rect from screen 0,0 to screen 1,1 fs solid noborder fc rgb "black"  
set label 5 at screen 0.5,0.97 center "{/Times-Italic=25 Atomic f-orbitals}" tc rgb "white"

set xrange [-0.75:0.75]
set yrange [-0.75:0.75]
set zrange [-0.75:0.75]
set pm3d depthorder
set view 60,140
set size ratio 1

unset key
unset colorbox
set border 0
unset xtics
unset ytics
unset ztics

hsize=0.21
aratio=10.0/7.0
vsize=hsize/aratio
hmarg=0.25
vmarg=0.08
vlabelmarg=0.02

bmarg=0.05
lmarg=(1.0-hmarg-2.0*hsize)*0.5
lmarg2 =(1.0-hsize)*0.5

set bmargin screen bmarg+(vsize+vmarg)*0
set tmargin screen bmarg+vsize*1+vmarg*0
set lmargin screen lmarg2
set rmargin screen lmarg2+hsize

Y(u,v)=Y30(u,v)
orbital="f_{z(5z^2-3r^2)} (f_{z^3})"
index = "l = 3, m = 0"
set label 1  at screen lmarg2+hsize*0.5,bmarg+vsize*0+vmarg*0-0.02 center sprintf("{/Time-Italic=18 %s}", orbital) tc rgb "white"
set label 10 at screen 0.5,bmarg+vsize*1+vmarg*0+vlabelmarg center sprintf("{/Time-Italic=20 %s}", index) tc rgb "white"

splot tablefile using (Fx($1,$2)):(Fy($1,$2)):(Fz($1,$2)):(light($1,$2)) with pm3d 

unset object
unset label 5
unset label 10

set bmargin screen bmarg+(vsize+vmarg)*1
set tmargin screen bmarg+vsize*2+vmarg*1
set lmargin screen lmarg+(hsize+hmarg)*0
set rmargin screen lmarg+hsize*1+hmarg*0

Y(u,v)=Y31_re(u,v)
orbital="f_{x(5z^2-r^2)} (f_{xz^2})"
index = "l = 3, m = {\261}1"
set label 1 at screen lmarg+hsize*0.5,bmarg+vsize*1+vmarg*1 center sprintf("{/Time-Italic=18 %s}", orbital) tc rgb "white"
set label 10 at screen 0.5,bmarg+vsize*2+vmarg*1+vlabelmarg center sprintf("{/Time-Italic=20 %s}", index) tc rgb "white"

splot tablefile using (Fx($1,$2)):(Fy($1,$2)):(Fz($1,$2)):(light($1,$2)) with pm3d 

unset label 10


set bmargin screen bmarg+(vsize+vmarg)*1
set tmargin screen bmarg+vsize*2+vmarg*1
set lmargin screen lmarg+(hsize+hmarg)*1
set rmargin screen lmarg+hsize*2+hmarg*1

Y(u,v)=Y31_im(u,v)
orbital="f_{y(5z^2-r^2)} (f_{yz^2})"
set label 1 at screen lmarg+hsize*1.5+hmarg,bmarg+vsize*1+vmarg*1 center sprintf("{/Time-Italic=18 %s}", orbital) tc rgb "white"

splot tablefile using (Fx($1,$2)):(Fy($1,$2)):(Fz($1,$2)):(light($1,$2)) with pm3d 


set bmargin screen bmarg+(vsize+vmarg)*2
set tmargin screen bmarg+vsize*3+vmarg*2
set lmargin screen lmarg+(hsize+hmarg)*0
set rmargin screen lmarg+hsize*1+hmarg*0


Y(u,v)=Y32_re(u,v)
orbital="f_{xyz}"
index = "l = 3, m = {\261}2"
set label 1 at screen lmarg+hsize*0.5+hmarg*0,bmarg+vsize*2+vmarg*2 center sprintf("{/Time-Italic=18 %s}", orbital) tc rgb "white"
set label 10 at screen 0.5,bmarg+vsize*3+vmarg*2+vlabelmarg center sprintf("{/Time-Italic=20 %s}", index) tc rgb "white"
splot tablefile using (Fx($1,$2)):(Fy($1,$2)):(Fz($1,$2)):(light($1,$2)) with pm3d 


unset label 10

set bmargin screen bmarg+(vsize+vmarg)*2
set tmargin screen bmarg+vsize*3+vmarg*2
set lmargin screen lmarg+(hsize+hmarg)*1
set rmargin screen lmarg+hsize*2+hmarg*1

Y(u,v)=Y32_im(u,v)
orbital="f_{z(x^2-y^2)}"

set label 1 at screen lmarg+hsize*1.5+hmarg,bmarg+vsize*2+vmarg*2 center sprintf("{/Time-Italic=18 %s}", orbital) tc rgb "white"
splot tablefile using (Fx($1,$2)):(Fy($1,$2)):(Fz($1,$2)):(light($1,$2)) with pm3d 

set bmargin screen bmarg+(vsize+vmarg)*3
set tmargin screen bmarg+vsize*4+vmarg*3
set lmargin screen lmarg+(hsize+hmarg)*0
set rmargin screen lmarg+hsize*1+hmarg*0

Y(u,v)=Y33_re(u,v)
orbital="f_{x(x^2-3y^2)}"
index = "l = 3, m = {\261}3"

set label 1 at screen lmarg+hsize*0.5+hmarg*0,bmarg+vsize*3+vmarg*3 center sprintf("{/Time-Italic=18 %s}", orbital) tc rgb "white"
set label 10 at screen 0.5,bmarg+vsize*4+vmarg*3+vlabelmarg center sprintf("{/Time-Italic=20 %s}", index) tc rgb "white"


splot tablefile using (Fx($1,$2)):(Fy($1,$2)):(Fz($1,$2)):(light($1,$2)) with pm3d 

unset label 10


set bmargin screen bmarg+(vsize+vmarg)*3
set tmargin screen bmarg+vsize*4+vmarg*3
set lmargin screen lmarg+(hsize+hmarg)*1
set rmargin screen lmarg+hsize*2+hmarg*1

Y(u,v)=Y33_im(u,v)
orbital="f_{y(3x^2-y^2)}"
set label 1 at screen lmarg+hsize*1.5+hmarg,bmarg+vsize*3+vmarg*3 center sprintf("{/Time-Italic=18 %s}", orbital) tc rgb "white"
splot tablefile using (Fx($1,$2)):(Fy($1,$2)):(Fz($1,$2)):(light($1,$2)) with pm3d 



unset multiplot

set size 1,1

N=4
print n
if (n==0) pause -1;\
     set out "atomic-f-orbitals.ps";\
     set term post portrait color solid size 7in,10in enhanced;\
     
if (n < N) ;\
     TERM="PS";\
     n=n+1;\
     reread

n=-1
TERM="DEFAULT"
set term pop
set out