# one way to plot 3D hydrogen atom orbitals

I wrote a code to plot 3D hydrogen atom orbitals. I first expressed the orbital wavefunction in spherical coordinates and then transformed to xyz coordinates in order to do isosurface Gizmo plot in Igor.

function create()
variable n=5,l=3,m=0
variable size1=400,size2=400
wave w=calcparametric(N,L,M,size1,size2)
execute "Gizmo0()"
end

function/wave calcparametric(N,L,M,size1,size2)
variable N,L,M,size1,size2
make/d/o/n=((size1)^2*(size2)) wavex,wavey,wavez,ff
variable i,j,k,dt,df,psi2,rr,theta,phi,t=0
variable xx,yy,zz
dt=pi/size1
df=2*dt
for(i=0;i<size1;i+=1)
theta=i*dt
for(j=0;j<size1;j+=1)
phi=j*df
for(k=0;k<size2;k+=1)
rr=k*0.05
variable temp1=factorial(n+l),temp2=factorial(n-l-1),temp3=laguerrea(n-l-1,2*l+1,2*rr/n/0.529)
psi2=magsqr(sqrt((2/n/0.529)^3/2/n/(temp1)^3*temp2)*exp(-rr/n/0.529)*(2*rr/n/0.529)^l*temp3*sphericalHarmonics(l,m,theta,phi))
ff[t]=psi2
wavex[t]=rr*sin(theta)*cos(phi)
wavey[t]=rr*sin(theta)*sin(phi)
wavez[t]=rr*cos(theta)
t+=1
endfor
endfor
endfor
concatenate/np=1 {wavex,wavey,wavez,ff},dest

Variable pee, qew, kappa
variable u,v,z
make/n=((size2)/10,(size2)/10,(size2)/10)/d/o kcube=0
setscale/p x,-0.05*size2,1,kcube
setscale/p y,-0.05*size2,1,kcube
setscale/p z,-0.05*size2,1,kcube
duplicate/o kcube, BlueJacobian
BlueJacobian = 0
for (j=0; j<DimSize(dest,0); j+=1)
pee = round((dest[j][0] - DimOffset(kCube,0))/DimDelta(kCube,0))
qew = round((dest[j][1] - DimOffset(kCube,1))/DimDelta(kCube,1))
kappa=round((dest[j][2] - DimOffset(kCube,2))/DimDelta(kCube,2))
if (pee >= 0 && pee < DimSize(kCube,0) && qew >= 0 && qew < DimSize(kCube,1)&& kappa >= 0 && kappa < DimSize(kCube,2))
kCube[pee][qew][kappa] += dest[j][3]
BlueJacobian[pee][qew][kappa] += 1
endif
endfor
for(u=0;u<DimSize(kCube,0);u+=1)
for(v=0;v<DimSize(kCube,1);v+=1)
for(z=0;z<dimsize(kcube,2);z+=1)
kCube[u][v][z] /= max(1,BlueJacobian[u][v][z])
endfor
endfor
endfor
killwaves dest,BlueJacobian,wavex,wavey,wavez,ff
return kcube
end

Window Gizmo0() : GizmoPlot
PauseUpdate; Silent 1       // building window...
// Building Gizmo 8 window...
NewGizmo/W=(21,44,330,320)
ModifyGizmo startRecMacro=700
ModifyGizmo scalingOption=63
AppendToGizmo isoSurface=root:kcube,name=isoSurface0
ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ surfaceColorType,1}
ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineColorType,0}
ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidthType,0}
ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ fillMode,2}
ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidth,1}
ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ isoValue,1e-08}
ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ frontColor,1,0,0,1}
ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ backColor,0,0,1,1}
ModifyGizmo setDisplayList=0, object=isoSurface0
ModifyGizmo autoscaling=1
ModifyGizmo currentGroupObject=""
ModifyGizmo showInfo
ModifyGizmo infoWindow={559,7,1392,321}
ModifyGizmo endRecMacro
ModifyGizmo SETQUATERNION={-0.681597,0.118549,-0.123937,0.711347}
EndMacro

A few comments in no particular order:

1.  I assume that you are familiar with the old example under File Menu->Example Experiments->Visualization->SphericalHarmonicsDemo.

2.  Both the code above and the demo experiment could benefit from multithreading as it took a long time to compute.

3.  Your isoValue choice of 1e-8 is curious.  To determine a better value I recommend executing:

Make/d/n=40 ddd=kCube[20][0][p]
display ddd

This should indicate that a more appropriate isoValue may be ~2e-15.