# hydrogen2sb.mac, corrected, from Steeb & Hardy Ch.26

### This script allows to find energy levels of Hydrogen atom from Schroedinger equation,
and here we illustrate how to use physical constants in Maxima

load (physical_constants);

Let us use notation from S.Fluegge "Practical Quantum Mechanic", vol.1, problem 61.
Kepler's problem for infinetely heavy point-like nucleus of charge Z. For Z=1 we have hydrogen atom.
Denote:

2*m0*E/hb^2 = -%gamma^2;

 here hb denotes \hbar, i.e. Planck constant divided by 2\pi.  Let

kappa = Z*e0^2*m0/(%gamma*hb^2);

Now for usual substitution in wavefunction u(r,theta,phi) to separate variables:

u = (1/r) * chi(r) * Ylm(theta,phi);

we write the Schroedinger equation as (with y in place of chi):

'diff(y, r, 2)  + (-%gamma^2 + 2*%gamma*kappa/r-l*(l+1)/r^2)* y = 0;

Define orbital number l, magnetic m, and radial k, e.g.:

l: 1; m: 0; k: 1;

now the principal quantum number n is

n:k+l+1;

alpha: e0^2/(hb*c);

alpha is fine structure constant.
Now substitute a polynomial v1 and exp as an ansatz for y, where

v1: (sum (a [j] *r^j, j, 1 ,k));

v1: v1+1;

y1: r^(l+1)*exp(-sqrt(g2)*r)*v1;

 y1 is the trial function of the radial wave function.

 Insert the ansatz into Schroedinger differential equation

left: diff(y1,r,2)+(-g2+2*m0*c*alpha/(hb*r)-l*(l+l)/(r*r))*y1;
left: left*exp(sqrt(g2)*r);
left: left/r;
left: ratsimp(left);
left: num(left);
left: expand(left);
left: subst(x,sqrt(g2),left);

h[0]: coeff(left,r,0); 
h[1]: coeff(left,r,1);
h[2]: coeff(left,r,2); 
h[3]: coeff(left,r,3);
h[4]: coeff(left,r,4);

s: solve ([h[0]=0, h[1]=0, h[2]=0, h[3] =0, h[4]=0],
   [x,a[1],a[2],a[3],a[4]]);

sol: last(s);

sol: first(sol);

sol: rhs(sol);

E_kl: -hb^2*sol*sol/2/m0;

 energy eigenvalue for k, l 

print ("E_kl=" ,E_kl);

The output is the energy eigenvalue E_kl.
Now we may use maxima package "physical_constant" which has a set of many dimensional constants.
It uses package ezunits which by default employs SI. 
One can try
demo(ezunits);
display_known_unit_conversions;

We can list the constants:

propvars (physical_constant);

And print some of them, e.g.

constvalue(%alpha);

constvalue (%c);

we may copy paste this SI value into CGS value for speed of light (adding two orders):

c0:29979245800;

constvalue (%h_bar);

float(%);

we may copy paste this SI value into CGS value (7 orders difference in J and erg):

hb:1.054571628251774e-27;

constvalue(%%e);

float(%);

e0:1.602176487e-19*c0/10;

constvalue(%m_e);

In the same way

m0:9.10938215e-28;

Now we get the energy of the level in ergs 

E: ev(E_kl);

and transform to electron-Volts (do not confuse evaluation operator ev above with eV)

load("physconst.mac");

float(%Ry);

Here we get Rydberg in inverse meters. Next line gives its uncertainty:

 get(%Ry, RSU);

We better get Ry in eV. We must multiply electron charge e0 by one Volt, expressed in CGS. The unit of voltage in CGS is 
1 statvolt = 299.792458 volts. 
(The conversion factor 299.792458 is simply the numerical value of the speed of light in cm/s divided by 10^8). Hence,

eV0:e0/(c0*1e-8);

This is the value of eV in ergs, now we have the energy of our level in eV:

E_eV:E/eV0;

Since Rydberg is

Ry0:m0*e0^4/hb^2;

in ergs, and in eV

Ry0eV:Ry0/eV0;

We obtain the energy of our level in Rydbergs:

E/Ry0;

It is easy to check that we get -1/(2n^2) when the principal quantum number of the level is n