BeginPackage["CentralField`"];
\(\(\(Clear[Radial, Angle, AnglePlot, Orbital];\)\(
\)
    \)\)
\(\(\(Radial::usage = "Radial[ro_, n_, l_, Z_] calculates the radial representation
of the eigenfunctions for an electron in the Coulomb potential. The numbers
n and l are the quantum numbers for the energy and the angular momentum
operator. Z specifies the number of charges in the nucleus. The radial
distance between the center and the electron is given by ro.";\)\(
\)
    \)\)
\(\(\(Angle::usage = "Angle[theta_, phi_, l_, m_] calculates the angular part of
the wave function for an electron in the Coulomb potential. The numbers L
and m denote the quantum numbers for the angular momentum operator. Theta
and phi are the angles in the spherical coordinate system.";\)\(
\)
    \)\)
\(\(\(Orbital::usage = "Orbital[theta_,phi_,l_,m_,type_String]  calculates the
superposition of two wave functions for the quantum numbers m_l = +m and
m_l = -m. The variable type allows the creation of the sum or the difference
of the wave functions. The string values of type are either plus or minus.";\)\(
\)
    \)\)
\(\(AnglePlot::usage = "AnglePlot[pl_,theta_,phi_] gives a graphical
representation of the function contained in pl. The range of representation
is Pi <= phi < 5 Pi/2 and 0 < theta < Pi. Theta is measured with
respect to the vertical axis. This function is useful for ploting the orbitals
and the angular part of the eigenfunction.";\)

 (*\(\(--\(-define\)\)\ global\ \(variables--\)\)\(-\)*) 
\)
theta::usage;
phi::usage;
ro::usage;
n::usage;
l::usage;
\(\(\(m::usage;\)\(
\)
    \)\)
\(\(Begin["`Private`"];\)

 (*\(\(--\(-radial\)\)\ part\ of\ the\ eigenfunctions\ in\ the\ Coulomb\ \(potential--\)\)\(-\)*) 
\)
\(\(Radial[ro_, n_, l_, Z_] := 
        Block[{norm, 
            hnl}, 
 (*\(--\(-\(normalization--\)\)\)\(-\)*) 
norm = \((Sqrt[\(\((n + 
                            l)\)!\)/\((2\ n\ \(\((n - l - 
                                1)\)!\))\)]\ \((\((2\ Z)\)/n)\)^\((l + 
                        3/2)\))\)/\(\((2\ l + 
                    1)\)!\); 
 (*\(\(--\(-definition\)\)\ of\ the\ wave\ \(function--\)\)\(-\)*) 
hnl = 
            norm\ ro^l\ Exp[\(-\((\((Z\ ro)\)/n)\)\)]\ Hypergeometric1F1[
                l + 1 - n, 
                2\ l + 2, \((2\ Z\ ro)\)/
                  n]];\)


 (*\(\(--\(-angular\)\)\ part\ of\ the\ eigenfunctions\ in\ the\ Coulomb\ \(field--\)\)\(-\)*) 
\)
Angle[theta_, phi_, l_, m_] := 
        Block[{norm, legendre, x, angle, m1, result}, 
          m1 = Abs[
              m]; 
 (*---normalization---*) 
norm = (-1)^
                m1 Sqrt[((2 l + 
                            1) (l - m1)!)/(2 (l + 
                              m1)!)]/
                Sqrt[2 Pi]; 
 (*---eigenfunctions---*) 
legendre = 
            Sin[theta]^m1 D[
                LegendreP[l, x], {x, m1}]; 
legendre = 
            legendre /. 
              x -> 
                Cos[theta]; 
 (*---consider the cases m > 0 and m < 0---*) 
If[
            m ≥ 0, angle = Exp[I m phi], 
            angle = (-1)^
                  m1 Exp[-(I m1 phi)]]; 
 (*---normalized eigenfunction---*) 
result = 
            norm legendre angle];

 (*--create orbitals---*) 

\(\(Orbital[theta_, phi_, l_, m_, type_String] := 
        Block[{norm, m1, rule, wave, wave2}, 
m1 = 
            Abs[m]; 
 (*\(\(--\(-replacement\)\)\ rule\ for\ the\ exponential\ \(function--\)\)\(-\)*) 
rule = {E^\((Complex[0, a_]\ \((x_. )\))\) -> 
                Cos[x\ Abs[a]] + 
                  I\ Sign[a]\ Sin[
                      x\ Abs[a]]}; 
 (*\(\(--\(-distinguish\)\)\ different\ \(cases--\)\)\(-\)*) 
If[
            m1 ≥ 1, 
If[
              type == "plus", 
 (*\(\(--\(-sum\)\)\ of\ the\ wave\ functions\ for\ a\ fixed\ \(m--\)\)\(-\)*) 
wave = 
                Expand[Angle[theta, phi, l, m1] + 
                      Angle[theta, phi, l, \(-m1\)] /. 
                    rule],  (*\(\(--\(-difference\)\)\ of\ the\ wave\ function\ for\ a\ fixed\ \(m--\)\)\(-\)*) 
wave = 
                Expand[Angle[theta, phi, l, m1] - 
                      Angle[theta, phi, l, \(-m1\)] /. 
                    rule]]; 
wave2 = 
              wave^2; 
 (*\(\(--\(-normalization\)\)\ of\ the\ \(superposition--\)\)\(-\)*) 
norm = 
              Integrate[
                wave2, {phi, 0, 2  Pi}, {theta, 0, 
                  Pi}]; 
wave2 = Expand[wave2/Abs[norm]], 
            wave = Angle[theta, phi, l, 
                  m1]^2]];\)

 (*\(\(--\(-graphical\)\)\ representation\ of\ the\ angular\ \(part--\)\)\(-\)*) 
\)
AnglePlot[pl_, theta_, phi_] := 
        Block[{}, 
  
            ParametricPlot3D[{-pl Sin[theta] Cos[phi], -pl Sin[
                  theta] Sin[phi], pl Cos[theta]}, {phi, Pi, 
              5  Pi/2}, {theta, 0, Pi}, PlotRange -> All, 
            PlotPoints -> {40, 40}]];
End[];
EndPackage[];

?AnglePlot

?Orbital

AnglePlot[
      Orbital[θ, ϕ, 4, 
        2, "minus"], θ, ϕ]

AnglePlot[
      Orbital[θ, ϕ, 4, 
        2, "minus"], θ, ϕ]