```In an [presentation][1] by Markus van Almsick, he gives an solution to visualize atomic orbitals using Image3D.

R[n_Integer?Positive, l_Integer?NonNegative, r_] :=
Block[{ρ = (2 r)/n},
Sqrt[(2/n)^3 (n - l - 1)!/(2 n (n + l)!)] E^(-ρ/2) ρ^l LaguerreL[n - l - 1, 2 l + 1, ρ]] /; l < n

full wave function:

ψ[n_, l_, m_, r_, ϑ_, φ_] := ψ[n, l, m, r, ϑ, φ] =
FullSimplify[R[n, l, r] SphericalHarmonicY[l, m, ϑ, φ], {r >= 0, ϑ ∈ Reals, φ ∈ Reals}]

CompileWaveFunction = Compile[{{x, _Real}, {y, _Real}, {z, _Real}},
Block[{ρ = x^2 + y^2, r, ϑ, φ},
If[ρ > 0,
r = Sqrt[ρ + z^2]; ϑ = ArcCos[z/r]; φ = ArcTan[x, y],
r = Abs[z]; ϑ = π/2 Sign[z]; φ = 0];
#
],
CompilationTarget -> "C"
] &;

color function:

colorFunction = (Blend[{
{0., RGBColor[0.7, 0.8, 1., 0.]},
{0.1, RGBColor[0., 0.7, 0.1, 0.012]},
{0.4, RGBColor[1., 0.1, 0.03169, 0.06723]},
{1., RGBColor[1., 0.95051, 0., 0.10963]}}, #] &)

plot 3p orbital

Block[
{nψ =
CompileWaveFunction[ψ[3, 1, 0, r, ϑ, φ]], data, vol},
data = Table[Abs[nψ[x, y, z]]^2, {z, -20, 20, 0.25}, {y, -20, 20, 0.25}, {x, -20, 20, 0.25}];
vol = RawArray["Byte", Round[(255/Max[data]) data]];
Image3D[vol, "Byte", Background -> Black,
Method -> {"FastRendering" -> True, "InterpolateValues" -> True},
ColorFunction -> colorFunction, BoxRatios -> 1]
]

![enter image description here][2]

it's quite smooth to rotate the Image3D object

![enter image description here][3]

----

We can also visualize the atomic orbital by plotting the isosurface:

Block[{nψ = CompileWaveFunction[ψ[3, 2, 0, r, ϑ, φ]]},
ContourPlot3D[Abs[nψ[x, y, z]]^2, {x, -20, 20}, {y, -20, 20}, {z, -20, 20},
PlotPoints -> 15, Contours -> {0.00002},
ColorFunction -> Function[{x, y, z}, ColorData["Rainbow"]@Rescale[Arg[nψ[x, y, z]], {-π, π}]],
ColorFunctionScaling -> False, Mesh -> None]
]

![enter image description here][4]

and make a animation showing different isosurfaces

plots = ParallelTable[
Block[{nψ = CompileWaveFunction[ψ[4, 2, 1, r, ϑ, φ]]},
ContourPlot3D[
Abs[nψ[x, y, z]]^2, {x, -20, 20}, {y, -20, 20}, {z, -20, 20}, PlotPoints -> 17, Contours -> {ct},
ColorFunction -> Function[{x, y, z}, ColorData["Rainbow"]@Rescale[Arg[nψ[x, y, z]], {-π, π}]],
Boxed -> False, Axes -> False, ColorFunctionScaling -> False, Mesh -> None,
ViewPoint -> {0.98, -2.76, 1.7}, ViewVertical -> {-0.004, -0.117, 0.993}]],
{ct, 0.00003, 0.000015, -0.0000005}
];
ListAnimate[plots]

![enter image description here][5]