Hydrogen Probability Density Movies

David McIntyre
13 March 2012

SetOptions[Plot,BaseStyle->{FontFamily->"Arial",FontSize->12},PlotStyle->{{Black,AbsoluteThickness[1.5]}},AxesStyle->Directive[FontFamily->"Arial",FontSize->12,Black,AbsoluteThickness[0.5],Arrowheads[0.04]],TicksStyle->{{Black,AbsoluteThickness[0.5]},{Black,AbsoluteThickness[0.5]}}];

Define Hydrogen Atom Eigenstate Wave Functions

(Subscript[R, n_,l_])[r_]:=(2/((n^2)(2l+1)!))(Sqrt[(((n+l)!)/((n-l-1)!))])Exp[-r/n]((((2/n)r))^l)Hypergeometric1F1[-n+l+1,2l+2,(2/n)r]
(Subscript[Ψ, n_,l_,m_])[r_,θ_,φ_]:=(Subscript[R, n,l])[r]SphericalHarmonicY[l,m,θ,φ]

Make Color Scale

colorscale=DensityPlot[y,{x,0,1},{y,0,1},PlotPoints->100,Mesh->False,AspectRatio->8,ImageSize->100,Axes->False,ColorFunction->(Hue[0.667 (1-#)]&),Frame->False,PlotRange->{{0,1},{-0.1,1.1}},Epilog->{Text[Style["0",12,FontFamily->"Arial"],{0.5,-0.05}],Text[Style["max",12,FontFamily->"Arial"],{0.5,1.05}]}]

Initialize Parameters

Execute this to choose nlm quantum numbers and to set up plot label, plot scale and number of plotting points (100 pts looks better, but is too slow for manipulate.  2 works OK)
Useful radial plot ranges: 1s: 3, 2s,2p:10, 3s,p,d:16

n=2;
l=1;
m=0;
lbl="2p"
rng=10;
Npts=20;
f[r_,θ_,ϕ_]:=(Subscript[Ψ, n,l,m])[r,θ,ϕ];
psimax=Max[Max[Table[(Abs[f[r,θ,0]]),{r,0.01,rng,0.1},{θ,0,π,π/20}]],Max[Table[(Abs[f[r,θ,π/2]]),{r,0.01,rng,0.1},{θ,0,π,π/20}]],Max[Table[(Abs[f[r,π/2,ϕ]]),{r,0.01,rng,0.1},{ϕ,0,2*π,π/40}]]]

Plot angular part

fang[θ_,ϕ_]:=SphericalHarmonicY[l,m,θ,ϕ];

Use color-on-sphere plot.  Red is max, blue is min (0).  Viewpoint is from positive x-axis.

maxx=Max[Max[Table[((Abs[N[fang[θ,0]]])^2),{θ,0,π,π/100}]],Max[Table[((Abs[N[fang[θ,π/2]]])^2),{θ,0,π,π/100}]],
Max[Table[((Abs[N[fang[π/2,ϕ]]])^2),{ϕ,0,2*π,π/100}]]]

ParametricPlot3D[{Cos[v] Sin[u],Sin[v] Sin[u],Cos[u]},{u,0,π},{v,0,2*π},PlotPoints -> 50,ColorFunction->Function[{x,y,z,u,v},Hue[0.67*(1-((Abs[fang[u,v]])^2)/maxx)]],ColorFunctionScaling->False,ViewPoint->{2.6, 0, 1.6},Axes->False,Boxed->False,Mesh->False,ImageSize->300]

Plot radial part

fr[r_]:=(Subscript[R, n,l])[r]

Plot[(((fr[r]))^2),{r,0,rng},AxesLabel->{"r", "R(r)^2"}]

Make static y = 0 slice

y=0.0;
yplot=ParametricPlot3D[{x,y,z},{x,-rng,rng},{z,-rng,rng},PlotPoints->Npts,PlotRange->{{-rng,rng},{-rng,rng},{-rng,rng}},AspectRatio->1,AxesLabel->{"x/(Subscript[a, 0])","y/(Subscript[a, 0])","z/(Subscript[a, 0])"},Mesh->None,ViewPoint->{5,3,2},AxesEdge->{{1,-1},{1,-1},{1,-1}},BaseStyle->{FontSize->12,FontFamily->"Arial"},ColorFunction->Function[{x,y,z},Hue[0.67* (1-((Abs[f[(Sqrt[(x^2)+(y^2)+(z^2)]),ArcCos[(z/(Sqrt[(x^2)+(y^2)+(z^2)]))],ArcCos[(x/(Sqrt[(x^2)+(y^2)]))]]])/psimax))]],ColorFunctionScaling->False]

Manipulate (animate) frames of z slices

Manipulate[
Graphics[{Inset[Show[ParametricPlot3D[{x,y,z},{x,-rng,rng},{y,-rng,rng},PlotPoints->Npts,PlotRange->{{-rng,rng},{-rng,rng},{-rng,rng}},AspectRatio->1,AxesLabel->{"x/(Subscript[a, 0])","y/(Subscript[a, 0])","z/(Subscript[a, 0])"},Mesh->None,ViewPoint->{5,3,2},AxesEdge->{{1,-1},{1,-1},{1,-1}},BaseStyle->{FontSize->12,FontFamily->"Arial"},ColorFunction->Function[{x,y,z},Hue[0.67* (1-((Abs[f[(Sqrt[(x^2)+(y^2)+(z^2)]),ArcCos[(z/(Sqrt[(x^2)+(y^2)+(z^2)]))],ArcCos[(x/(Sqrt[(x^2)+(y^2)]))]]])/psimax))]],ColorFunctionScaling->False],yplot],{0,0},Center,{2,2}],Inset[colorscale,{1.1,0.25},Center,{1,1}],Inset[Style[lbl,14,FontFamily->"Arial"],{-0.7,0.86},Center,Automatic]},ImageSize->488],{z,-rng,rng,rng/20.0}]