#
#   Algebraic Function Approximant Program
#
#	by G.J. Fee and 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 a specified state where n <> |k|
#
Edirac:=mc2/sqrt( (1 + (Z*alp/(n-k + k*sqrt(1-(Z*alp/k)^2)))^2) );
# 
#   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;
#
#   Try "4s" state
#
n:=4; k:=1;
#
Etrial := Edirac;
#
#   Assume only even powers of alp and y .
#
Aeq := c00 + c20*alp^2 + c02*y^2 + c40*alp^4 + c22*alp^2*y^2 + c04*y^4 
       + c42*alp^4*y^2 + c24*alp^2*y^4;
Order := 16;
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)/c40);
quit;
