Update:
Simpler version of the code with ColorFunctionScaling -> {0, Pi^2}
and for some reason also plot points had to be increased to PlotPoints -> 300
to achieve the same quality.
NormFactor[n_, l_] := Sqrt[Factorial[n - l - 1]/(2 n Factorial[n + l])]
Psi[n_, l_, m_, u_, \[Theta]_, \[Phi]_] :=
NormFactor[n, l] Exp[-u/2] u^l LaguerreL[n - l - 1, 2 l + 1,
u] SphericalHarmonicY[l, m, \[Theta], \[Phi]]
u[x_, z_] := Sqrt[x^2 + z^2]
theta[x_, z_] := ArcTan[z, x]
With[{nn = 3, nl = 1, nm = 0},
DensityPlot[
Abs[Psi[nn, nl, nm, u[x, z], theta[x, z], 0]]^2, {x, -15,
15}, {z, -15, 15}, Frame -> None, Axes -> False, PlotRange -> All,
Ticks -> None, PlotPoints -> 300, ColorFunctionScaling -> {0, Pi^2},
ColorFunction -> "SunsetColors", ImageSize -> {400, 400}]]
With ColorFunction -> GrayLevel
:
With ColorFunction -> (GrayLevel[1 - #] &)
:
Old version:
So I think that they are simply clipping the range of the color gradient (that also explain the big white areas on their images).
The constant 0.005637338004403787
is the max value of the function and the clipping constant is Pi^2
or could be just, say, 9
or 10
. Also you have to use ColorFunctionScaling -> False
for a proper effect of clipping of the gradient.
With[{nn = 3, nl = 1, nm = 0},
DensityPlot[Pi^2/
0.005637338004403787` Abs[
Psi[nn, nl, nm, u[x, z], theta[x, z], 0]]^2, {x, -15,
15}, {z, -15, 15}, Frame -> None, Axes -> False, PlotRange -> All,
Ticks -> None, PlotPoints -> 120, ColorFunctionScaling -> False,
ColorFunction -> "SunsetColors", ImageSize -> {400, 400}]]
My image is almost unrecognizable form the one on Wikipedia:
With ColorFunction -> ColorData[{"GrayTones", "Reverse"}]
you get: