(* This notebook attempts to make 3D plots of user selected hydrogen states *)

Clear[a,r,θ,ϕ,x,y,z,n,l,m,Lq,Lpq,ψnml]

(* Select Wavenumbers *)
n=5;
l=4; (* (Subscript[l, max])= n-1 *)
m=4;(* |m| ≤ l *)

(* Generate Radial Equation *)
p=2l+1;
q=p+n-l-1;
(* Leguerre Polynomials *)
Lq[x_]=Exp[x]D[Exp[-x]x^q,{x,q}]; 
Lpq[x_]=(-1)^p D[Lq[x],{x,p}];
Rnl[r_]=FullSimplify[Sqrt[(2/(n a))^3 (n-l-1)!/(2n((n+l)!)^3)]Exp[-r/(n a)](2r/(n a))^l Lpq[2r/(n a)]];

(* Generate Angular equation *)
Pl[x_]=Simplify[1/(2^l l!)D[(x^2-1)^l,{x,l}]];
Plm[x_]=(1-x^2)^(Abs[m]/2)D[Pl[x],{x,Abs[m]}];
ϵ=Piecewise[{{(-1)^m,m≥0},{1,m<0}}];
(* Legendre Polynomials *)
Ylm[x_,ϕ_]=ϵ Sqrt[((2l+1)(l-Abs[m])!)/(4 Pi (l+Abs[m])!)] Exp[I m ϕ]Plm[x]; 
Ylm[θ_,ϕ_]=FullSimplify[Ylm[Cos[θ],ϕ]];

(* Total Wave Equation *)
ψnml[r_,θ_,ϕ_]=FullSimplify[Rnl[r]Ylm[θ,ϕ],{Sin[θ]>0,Cos[θ]>0}]

(* a is the Bohr Radius. a = 4(Subscript[Piϵ, 0])𝒽/me^2 = 0.529 10^-10 m *)
a=1;
Clear[x,y,z]
(* Convert to Cartesian *)
Mess=Solve[{x==r Sin[θ]Cos[ϕ],y==r Sin[θ]Sin[ϕ],z==r Cos[θ]},{r,θ,ϕ}];
R=r/.Mess[[16]];
Θ=θ/.Mess[[16]];
Φ=ϕ/.Mess[[16]];

Density[r_,θ_,ϕ_]=FullSimplify[Re[ψnml[r,θ,ϕ]^2],{r∈Reals,θ∈Reals,ϕ∈Reals}];
DensityCart[x_,y_,z_]=FullSimplify[Density[R,Θ,Φ]]

S=45;
Q=.85;
Thresh=0.000003;

DensityPlot[DensityCart[x,0,z],{x,-S,S},{z,-S,S},PlotPoints->150,ColorFunction->"GrayTones"]

RegionPlot3D[DensityCart[x,y,z]>Thresh,{x,-Q S,Q S},{y,-Q S,Q S},{z,-Q S,Q S},Boxed->False,Axes->False,PlotPoints->50,PlotRange->Full]

{n,l,m}

(* Saves density plot and animated rotation about z axis *)


axis={0,0,1};
s=Table[RegionPlot3D[DensityCart[x,y,z]>Thresh,{x,-Q S,Q S},{y,-Q S,Q S},{z,-Q S,Q S},Boxed->False,Axes->False,PlotPoints->50,PlotRange->Full,SphericalRegion->True]/.gg:GraphicsComplex[___]:>Rotate[gg,theta,axis],{theta,0.,2. Pi(19/20),2Pi/20}];

Export["544.gif",s]
Export["544.jpg",

DensityPlot[DensityCart[x,0,z],{x,-S,S},{z,-S,S},PlotPoints->150,ColorFunction->"GrayTones",Axes->False]
]