################################################################################
# make polar plots of real combinations of spherical harmonics                 #
# > gnuplot -persist orbitals.gnu                                              #
#                                                    I 2001, E. Koch,  MPI-FKF #
################################################################################

# set term postscript portrait enh 'Helvetica' 21
# set output 'tmp.ps'
# set size 1.0, 0.7

# set up coordinate system #####
set view 60, 30
scale=0.9
set xrange [-scale:scale]
set yrange [-scale:scale]
set zrange [-scale:scale]
set noborder
set noxtics
set noytics
set noztics
set nokey
# coordinate axes
set arrow 1 from -scale, 0, 0 to scale, 0, 0 lw 3 
set arrow 2 from 0, -scale, 0 to 0, scale, 0 lw 3
set arrow 3 from 0, 0, -scale to 0, 0, scale lw 3
set label 'x' at scale,-.25*scale,0
set label 'y' at 0,1.1*scale,-.15*scale
set label 'z' at 0.10*scale,0,scale

# polar coordinates ####
set parametric
set urange [0:pi]
set vrange [-pi:pi]
x(u,v)=sin(u)*cos(v)
y(u,v)=sin(u)*sin(v)
z(u,v)=cos(u)
set isosamples 22, 45
set hidden3d

# atomic orbitals ####
# s orbital
s(u,v)=sqrt(1.0/4.0/pi)
# p orbitals
pz(u,v)=sqrt(3.0/4.0/pi)*cos(u)         #  m =0
px(u,v)=sqrt(3.0/4.0/pi)*sin(u)*cos(v)  # |m|=1
py(u,v)=sqrt(3.0/4.0/pi)*sin(u)*sin(v)  # |m|=1
# d orbitals
d3z2m1(u,v)=sqrt( 5.0/16.0/pi)*(3.0*cos(u)**2-1.0)           #  m =0
dzx(u,v)   =sqrt(15.0/16.0/pi)*sin(2.0*u)*cos(v)             # |m|=1
dyz(u,v)   =sqrt(15.0/16.0/pi)*sin(2.0*u)*sin(v)             # |m|=1
dx2my2(u,v)=sqrt(15.0/16.0/pi)*sin(u)**2*cos(2.0*v)          # |m|=2
dxy(u,v)   =sqrt(15.0/16.0/pi)*sin(u)**2*sin(2.0*v)          # |m|=2
# f orbitals
fz5z2m3 (u,v)=sqrt(  7.0/16.0/pi)*(5.0*cos(u)**2-3.0)*cos(u)        #  m =0
fx5z2m1 (u,v)=sqrt( 21.0/32.0/pi)*(5.0*cos(u)**2-1.0)*sin(u)*cos(v) # |m|=1
fy5z2m1 (u,v)=sqrt( 21.0/32.0/pi)*(5.0*cos(u)**2-1.0)*sin(u)*sin(v) # |m|=1
fzx2my2 (u,v)=sqrt(105.0/16.0/pi)*cos(u)*sin(u)**2*cos(2.0*v)       # |m|=2
fxyz    (u,v)=sqrt(105.0/16.0/pi)*cos(u)*sin(u)**2*sin(2.0*v)       # |m|=2
fxx2m3y2(u,v)=sqrt( 35.0/32.0/pi)*sin(u)**3*cos(3.0*v)              # |m|=3
fy3x2my2(u,v)=sqrt( 35.0/32.0/pi)*sin(u)**3*sin(3.0*v)              # |m|=3

# some hybrid orbitals 
# sp hybrids:
sp_1(u,v)=sqrt(1.0/2.0)*(s(u,v)+pz(u,v))
sp_2(u,v)=sqrt(1.0/2.0)*(s(u,v)-pz(u,v))
# sp2 hybrids:
sp2_1(u,v)=sqrt(1.0/3.0)*(s(u,v)+sqrt(2.0)*px(u,v))
sp2_2(u,v)=sqrt(1.0/3.0)*(s(u,v)-sqrt(1.0/2.0)*px(u,v)+sqrt(3.0/2.0)*py(u,v))
sp2_3(u,v)=sqrt(1.0/3.0)*(s(u,v)-sqrt(1.0/2.0)*px(u,v)-sqrt(3.0/2.0)*py(u,v))
# sp3 hybrids:
sp3_1(u,v)=0.5*(s(u,v)+px(u,v)+py(u,v)+pz(u,v))
sp3_2(u,v)=0.5*(s(u,v)+px(u,v)-py(u,v)-pz(u,v))
sp3_3(u,v)=0.5*(s(u,v)-px(u,v)+py(u,v)-pz(u,v))
sp3_4(u,v)=0.5*(s(u,v)-px(u,v)-py(u,v)+pz(u,v))

# >>> here you choose what orbital to plot <<< #################################
orb(u,v)=sp3_1(u,v) 
# plot
pos(u,v)= orb(u,v)<0 ? 0 :  orb(u,v)
neg(u,v)= orb(u,v)>0 ? 0 : -orb(u,v) 
splot pos(u,v)*x(u,v), pos(u,v)*y(u,v), pos(u,v)*z(u,v), \
      neg(u,v)*x(u,v), neg(u,v)*y(u,v), neg(u,v)*z(u,v)