(* :Title: HydrogenicAtom.m *)

(* :Context: HydrogenicAtom` *)

(* :Author: Michael A. Morrison *)

(* :Summary:
   A package to generate the spatial functions for stationary
   angular momentum eigenstates of the hydrogenic atom.
 *)

(* :Copyright: © 1997 by Michael A. Morrison *)

(* :Package Version: 3.4 *)

(* :Mathematica Version: 3.0 *)

(* :History:
   Version 1 written for Physics 4803, 1993;
   Version 2 used in RA+Rg research, 1993--1996;
   Version 3.0 & 3.1 for Physics 4803, October 1997;
   Version 3.2 for Physics 4803, October 1998;
   Version 3.3: added j, momRadFn[p], theta and phi fns. December 88
   Version 3.4: export n, l, and ml (August 1999)
*)

(* :Keywords: hydrogen atom, Laguerre polynomials, spherical
harmonics package *)

(* :Sources:
   The Joy of Quantum Physics: A Second User's Guide, Chap. 5
*)

(* :Warnings:
 At present commands for statistical properties of the hydrogenic atom
are not "idiot proof." The user MUST impose correct restrictions on n,
l, ml, and Z.*) The package redefines Conjugate to handle expressions.

(* :Limitations:
   <special cases not handled, known problems>
*)

(* :Discussion:
The conventions of Messiah for Laguerre polynomials are currently
used. Version 3.1 removes the restriction to atomic units that was in
force in Version 3.0. Version 3.2 improves documentation and
incorporates modified Conjugate, Re, and Im commands. Conventions of
Bethe and Salpeter are used for the momentum space radial functions.
*)

(* :Requirements:
  The Mathematica definition of Conjugate must be generalized,
as in Phys4803`QuickReIm`.*)

(* :Examples: *)

(* set up the package context, including public imports *)

BeginPackage["HydrogenicAtom`"]

(* usage messages for the exported functions and context itself *)

HydrogenicAtom::usage = "HydrogenicAtom.m is a package for calculating
normalized radial function R_{n,l}(r) for the hydrogenic atom.  This
package uses the conventions of Merzbacher, Messiah, etc. to conform
to Mathematica's definition of the associated Laguerre polynomials.
The first Bohr radius a0 appears in wave functions, energies, and
expectation values; to get rid of it, set a0->1. See also the command
atomicUnits. To get a list of commands in this package, type
?HydrogenicAtom`*."

uFnH::usage = "uFnH[n,l,Z][r] is the reduced radial function for the
hydrogenic atom with principal quantum number n and orbital angular
momentum quantum number l. The default value for Z=1, atomic
hydrogen."

radFnH::usage = "radFnH[n,l,Z][r] is the radial function  in atomic
units for the hydrogenic atom with principal quantum number n and
orbital angular momentum quantum number l. The default value for Z=1,
atomic hydrogen."

momRadFnH::usage = "momRadFnH[n,l][p] is the momentum space radial function
in atomic units for the hydrogenic atom with principal quantum number n and
orbital angular momentum quantum number l. The present encoding works only
for atomic hydrogen (Z=1) and does not explicitly show a0."

radFnHrho::usage = "radFnHrho[n,l,Z][rho] is the radial function  in
atomic units for the hydrogenic atom with principal quantum number n
and orbital angular momentum quantum number l expressed as a function
of the auxilliary variable rho=2Zr/(na0). The default value for Z=1,
atomic hydrogen."

Z::usage = "Z is the nuclear charge in atomic units. The default value
for Z=1, atomic hydrogen. If you don't want the default, feed the symbol Z
or some value as the final argument in the function being calculated."

a0::usage = "a0 is the first Bohr radius for atomic hydrogen."

atomicUnits::usage = "atomicUnits is a list of replacement rules that
converts an expression written in terms of hbar, me, e, a0, Eh, and/or
Ry to atomic units."

eH::usage = "eH[n,Z] is the hydrogenic atom energy in Hartree. The
default value for Z=1, atomic hydrogen."

meanrH::usage = "meanrH[n,l,Z] is the mean radius (in Bohr) in an
(n,l) stationary state. The default value for Z=1, atomic hydrogen."

meanrsqH::usage = "meanrsqH[n,l,Z] is the mean of the radius squared
(in square Bohr) in an (n,l) stationary state. The default value for
Z=1, atomic hydrogen."

uncrH::usage = "uncrH[n,l,Z] is the position uncertainty (standard deviation)
in Bohr) in an (n,l) stationary state. The default value for Z=1,
atomic hydrogen."

uncpH::usage = "uncpH[n,l,Z] is the momentum uncertainty (standard deviation)
in Bohr) in an (n,l) stationary state. The default value for Z=1,
atomic hydrogen."

Y::usage = "Y[l,ml][\[Theta],\[Phi]] is the spherical harmonic."

ThetaFn::usage = "ThetaFn[l,ml][\[Theta]] is the normalized \[Theta] function in
the spherical harmonic."

PhiFn::usage = "PhiFn[ml][\[Phi]] is the normalized \[Phi] function in
the spherical harmonic."

Y::usage = "Y[l,ml][theta,phi] is the spherical harmonic."

rctpH::usage = "rctpH[n, \[ScriptL], Z] generates the classical turning point
(in Bohr) for a SAME."

psiH::usage = "psiH[n,l,ml,Z][r,theta,phi] is the normalized
hydrogenic spatial function for principal quantum number n, orbital
angular momentum quantum number l, and azimuthal quantum number ml.
The default value for Z=1,  atomic hydrogen."

AbsYlmSq::usage = "AbsYlmSq[l,m][theta] produces the squared modulus
of the spherical harmonic Ylm, a function of theta only."

jH::usage = "jH[n, \[ScriptL], ml][r, \[Theta]] generates the probability
current density (in atomic units) for a stationary angular momentum eigenstate
of atomic hydrogen."

n::usage = "n is the principal quantum number."
l::usage = "l is the orbital angular momentum quantum number."
ml::usage = "ml is the magnetic (azimuthal) quantum number."

Begin["`Private`"]    (* begin private context (implementation) *)

protected = Unprotect[ Re, Im, Conjugate ]  (* from QuickReIm` *)

Conjugate[expr_] := expr /. Complex[x_,y_] -> x - I y
Re[expr_] := Expand[ (expr + Conjugate[expr])/2    ]
Im[expr_] := Expand[ (expr - Conjugate[expr])/(2I) ]

(* definition of auxiliary functions and local (static) variables *)

(* definition of the exported functions *)

atomicUnits = {hbar->1, me->1, e->1, a0->1, Eh->1, Ryd->1/2};

Y[l_,ml_][theta_,phi_] := SphericalHarmonicY[l,ml,theta,phi]

ThetaFn[\[ScriptL]_, m_][\[Theta]_] := Sqrt[2*Pi]*Y[\[ScriptL], m][\[Theta], 0]

PhiFn[m_][\[Phi]_] := 1/Sqrt[2*Pi]*Exp[I*m*\[Phi]]

eH[n_,Z_:1] := -1/(2n^2) Z^2;

radFnHrho[n_, l_,Z_:1][rho_] :=
2/n^2 (Z/a0)^(3/2) Sqrt[(n-l-1)!/(n+l)!] *
rho^l * Exp[-rho/2] *
LaguerreL[n-l-1, 2*l+1, rho]

radFnH[n_, l_, Z_:1][r_] :=
2/n^2 (Z/a0)^(3/2) Sqrt[(n-l-1)!/(n+l)!] *
((2/n)(Z r/a0))^l * Exp[-(Z r/(n a0))] *
LaguerreL[n-l-1, 2*l+1, (2 Z r/(n a0))]

uFnH[n_, l_, Z_:1][r_] := radFnH[n,l,Z][r]*r

psiH[n_,l_,ml_,Z_:1][r_,theta_,phi_] :=
   radFnH[n,l,Z][r] * Y[l, ml][theta,phi];

meanrH[n_, l_,Z_:1] :=
(3*n^2 - l*(l + 1))/(2Z) a0;

meanrsqH[n_, l_,Z_:1] :=
n^2/(2Z^2)(5n^2+1-3*l*(l+1)) a0^2;

rctpH[n_, \[ScriptL]_, Z_:1] :=
  {(n*(n - Sqrt[n^2 - \[ScriptL]*(\[ScriptL] + 1)]))/
    Z, (n*(n + Sqrt[n^2 - \[ScriptL]*(\[ScriptL] + 1)]))/Z} a0;

uncrH[n_, l_,Z_:1] :=
Sqrt[n^2(n^2+2)-l^2(l+1)^2]/(2Z) a0;

uncpH[n_, \[ScriptL]_,z_:1] :=
  Z/a0/n*Sqrt[1 - (2*\[ScriptL]*(\[ScriptL] + 1))/(n*(2*\[ScriptL] + 1))];

AbsYlmSq[l_, ml_][theta_] := Module[ {phi},
    Abs[ SphericalHarmonicY[l,ml,theta,phi] /. phi->0]^2 ];

jH[n_, \[ScriptL]_, ml_][r_, \[Theta]_] :=
   ml/(r*Sin[\[Theta]])*psiH[n, \[ScriptL], ml][r, \[Theta], 0] /. atomicUnits

momRadFnH[n_, \[ScriptL]_][p_] :=
    p*Sqrt[2/Pi*(n - \[ScriptL] - 1)!/(n + \[ScriptL])!]*
    n^2* 2^(2*(\[ScriptL] + 1))*\[ScriptL]!*
   (n^\[ScriptL]*p^\[ScriptL])/ (n^2*p^2 + 1)^(\[ScriptL] + 2)*
  (GegenbauerC[n - \[ScriptL] - 1, \[ScriptL] + 1, x] /.
       x -> (n^2*p^2 - 1)/(n^2*p^2 + 1));

End[ ]         (* end the private context *)

Protect[ Release[protected] ] (* protect exported symbols *)

EndPackage[ ]  (* end the package context *)