share/orthopoly/h_atom.dem

   1 if get('orthopoly,'version) = 'false then load("orthopoly")$
   2 showtime : all;
   3
   4 /* Let %a0 be the Bohr radius, mu the reduced mass, and %eo the
   5 electron charge. Declare these constants to be positive. */
   6
   7 assume(%a0 > 0, mu > 0, %hbar > 0, %eo > 0);
   8
   9 /* Merzbacher and orthopoly use different normalizations for the
  10 generalized Laguerre polynomials.  The function convert_to_mertz
  11 supplies a factor that converts from orthopoly normalization to
  12 Merzbacher. The hydrogen atom  eigenfunctions  are given by */
  13
  14 convert_to_mertz(n,l) := ((n+l)!^2 / ((n-l-1)! * (2*l+1)!)) / binomial(n+l, n-l-1);
  15
  16 psi(n,l,m) := sqrt((2 / (%a0 * n))^3 * (n-l-1)! / (2*n * ((n+l)!)^3)) *
  17      (2 * r / (%a0 * n))^l * exp(-r / (%a0 * n)) * gen_laguerre(n-l-1,2*l+1,2*r / (%a0 * n))  *
  18       spherical_harmonic(l,m, theta,phi) *  convert_to_mertz(n,l) ;
  19
  20
  21 /* Define the L2 inner product using a matchfix operator <<f,g>>. */
  22
  23 matchfix("<<",">>");
  24
  25 "<<"(f,g) := integrate(integrate(integrate(conjugate(f)*g*r^2 * sin(theta),
  26    theta,0,%pi),phi,0,2*%pi),r,0,inf);
  27
  28 /* Find <r> ,<r^2>, and <r^2> -<r>^2 for n=1,2,3, and 4 states. */
  29
  30 makelist(makelist(<<psi(n,l,0), r * psi(n,l,0)>>,l,0,n-1),n,1,4);
  31 makelist(makelist(<<psi(n,l,0), r^2 * psi(n,l,0)>>,l,0,n-1),n,1,4);
  32 makelist(makelist(<<psi(n,l,0), r^2 * psi(n,l,0)>> -
  33    <<psi(n,l,0), r * psi(n,l,0)>>^2,l,0,n-1),n,1,4);
  34
  35 /* Find the energies of the n=1,2,3, and 4 states.  Let ham
  36 be the hamiltonian. The function replace_hbar(e) replaces
  37 hbar by %eo * sqrt(mu * %a0).*/
  38
  39 ham(e) := (-%hbar^2 / (2*mu)) * (diff(r^2 * diff(e,r),r) / r^2 + diff(sin(theta)*diff(e,theta),theta)/(r^2 * sin(theta))
  40  + diff(e,phi,2)/(r^2 * sin(theta)^2)) - %eo^2 *e / r;
  41
  42 replace_hbar(e) := radcan(subst(%hbar=%eo * sqrt(mu * %a0),e));
  43
  44 replace_hbar(makelist(makelist(<<psi(n,l,0), ham(psi(n,l,0))>>,l,0,n-1),n,1,4));
  45
  46 /* Use degenerate perturbational methods to study the n = 3 Stark
  47 hamiltonian; this is problem 17.5 in Merzbacher; alpha is the
  48 reduced electric field strength. */
  49
  50 n : 3;
  51
  52 e : makelist(makelist(psi(n,l,m),m,-l,l),l,0,n-1)$
  53
  54 e : apply(append,e)$
  55
  56 stark_pot : alpha * r * cos(theta);
  57
  58 overlap(f,g) := <<f, stark_pot * g>>;
  59
  60 m : apply(matrix, outermap(overlap, e,e))$
  61
  62 load("eigen");
  63
  64 energy_shifts : eigenvalues(m);