#
#   Algebraic Function Approximant Program
#
#	by G.J. Fee, T.C. Scott
#
#	This routine exactly reconstucts the expression of the
#	Dirac energies from the series coefficients using a
#	heuristic polynomial fit.
#	
#	This routine handles the special case where n = |k|, 
#	e.g.: ground state ( n=1 k=1 )
#
#
Edirac:=mc2/sqrt( (1 + (Z*alp/(n-k + k*sqrt(1-(Z*alp/k)^2)))^2) );
n := k;
# 
#   To make life simpler absorb atomic charge Z into alp and set the
#   rest mass to 1.  Energy will be in units of rest energy.
#
mc2:=1; Z := 1;
#
Etrial := Edirac;
#
#   Assume only even powers of alp and y.
#
Aeq := c00 + c01*y + c20*alp^2 + c02*y^2;
#
Order := 5;
T := taylor(Etrial,alp);
subs(y=T,Aeq):
taylor(",alp):
convert(",polynom):
coeffs(",alp);
sol := solve({"});
lprint(`The Energy is one of the roots of this expression`);
expand(subs(sol,Aeq)/c20);
quit;
