(* Comparison with ground state of helium *)
(* Model V(r) = -1/r + Lambda[1-exp(-Delta r)]/r *)
(*SetDirectory["~/math/screened"];*)
o="sum.out";
elist=ReadList[o,{Number,Number,Number,Number,Number}];
Dimensions[elist];
npoints=Dimensions[elist][[1]];
delta=0.88;
redat=imdat={};
Do[{lambdan,deltan,re,im,acc}=elist[[n]];
   If[deltan!=delta,Goto[1]];
   If[acc<3,Return[]];
   redat=Append[redat,{lambdan,-re}];
   imdat=Append[imdat,{lambdan,im}];
   Label[1],{n,npoints}];
      (* 1/Z-expansion for the ground state of helium *)
(*SetDirectory["~/math/he1z"];*)
ndegree=4; (* degree of algebraic approximant *)
lamin=0;
lamax=13/10;
lastep=1/100;
nen=2; (* only the last nen approximants will be calculated *)
o="he1z.dat";
elist=ReadList[o,{Number,Number}];
Dimensions[elist];
nm=Dimensions[elist][[1]];
Do[en[n]=elist[[n,2]],{n,nm}];
              (* Summation 0 *)
en0=en[1];
ndegall={ndegree};
ndig=32;
ndegm=Length[ndegall];
      Clear[cf,ene];
e1=Table[Null,{n,2,nm}];
a0=en[1];
Do[e1[[n-1]]=en[n],{n,2,nm}];
  Clear[en];
km0=2;
kmm=ndegall[[ndegm]]+1;
ek=Table[Null,{k,1,kmm},{n,1,nm}];
c=Table[Null,{k,1,kmm},{n,1,nm}];
r=Table[Null,{k,1,kmm},{k,1,kmm}];
   Do[ndegree=ndegall[[ndegn]];
Print["--------- degree = ",ndegree];
km=ndegree+1;
ek[[1,1]]=1;
Do[ek[[1,n]]=0, {n,2,nm}];
Do[ek[[2,n]]=e1[[n-1]], {n,2,nm}];
If[km<3,Goto[1]];
Do[Print["Evaluating E^",k-1];
   Do[ek[[k,n]]=Sum[e1[[n-l]]*ek[[k-1,l]],{l,k-1,n-1}],
      {n,k,nm}],{k,km0,km}];
Label[1];km0=km+1;
Do[Do[r[[k,n]]=(-a0)^(n-k)*Binomial[n-1,k-1],
      {n,1,km}],{k,1,km}];
Do[Do[c[[k,n]]=ek[[k,n]],{n,k,nm}],{k,1,km}];
Do[nf=ka-km;
   nf2=nf+2;
   Do[n=nm-n5+nf2;
      c[[1,n]]=c[[1,n-1]],{n5,nf2,nm}];
   c[[1,nf+1]]=1;
   Do[n0=nf+k;
      a=c[[1,n0]]/c[[k,n0]];
      c[[k,nf+1]]=-a;
      cf[k,nf+1]=N[c[[k,nf+1]],ndig];
      Do[c[[1,n]]=c[[1,n]]-a*c[[k,n]], {n,n0,nm}], {k,2,km}];
   Do[a=c[[1,n]];
      Do[c[[k-1,n]]=c[[k,n]], {k,2,km}];
   c[[km,n]]=a, {n,nf2,nm}], {ka,km,nm}]
,{ndegn,1,ndegm}];
Clear[e1,ek,c,r];
              (* Summation *)
acc=1; (* accuracy limit *)
roots=Table[Null,{n,1,nen},{nd,1,ndegree}];
r=Table[Null,{k,0,ndegree},{k,0,ndegree}];
energy=Table[Null,{n,1,nen}];
delta=Table[Null,{n,1,nen}];
km=ndegree+1;
Clear[p,pa,pb,pc,pd,pe,t];
p={pa,pb,pc,pd,pe};
rs=Roots[Sum[p[[n+1]]*t^n,{n,0,ndegree}]==0,t,Cubics->True,Quartics->True];
Do[
Do[Do[r[[k,n]]=(-a0)^(n-k)*Binomial[n-1,k-1],
      {n,1,km}],{k,1,km}];
x=SetPrecision[rla,ndig];
Do[n=ka;
   na=n-km+1;
   cf[1,na]=1;
   Do[a=x*cf[1,na]*r[[k,1]];
      Do[a=a+cf[m,na]*r[[k,m]],{m,2,km}];
      Do[r[[k,m-1]]=r[[k,m]],{m,2,km}];
      r[[k,km]]=a, {k,1,km}];
   If[n<nm-nen+1,Goto[1]];
   pa=r[[1,km]];pb=r[[2,km]];pc=r[[3,km]];pd=r[[4,km]];pe=r[[5,km]];
   rs0=N[rs,ndig];
   s0=Solve[rs0,t];
   s=N[s0];
   Do[roots[[nm-n+1,k]]=Replace[t,s[[k]][[1]]],{k,1,ndegree}];
   Label[1], {ka,km,nm}];
d=10^66;
Do[Do[d1=Abs[roots[[1,k1]]-roots[[2,k2]]];
   If[d1>d,Goto[3]];
   d=d1;
   k10=k1;k20=k2;
   Label[3],{k1,1,ndegree}],{k2,1,ndegree}];
y0=roots[[1,k10]];
Do[d=10^66;
  Do[d1=Abs[roots[[n,k]]-y0];
     If[d1>d,Goto[4]];
     k0=k;
     d=d1;
     Label[4],{k,1,ndegree}];
   energy[[n]]=roots[[n,k0]];
   delta[[n]]=d,{n,1,nen}];
d=0;
Do[d1=delta[[n]];
   If[d1<d,Goto[5]];
   d=d1;
   Label[5],{n,1,nen}];
If[d>acc,Exit[]];
yre=SetPrecision[Re[y0],5];
yim0=Abs[SetPrecision[Im[y0],5]];
yim=0;
If[yim0>1.1*10^-5,yim=yim0];
d=SetPrecision[-Log[10,d],2];
ene[rla]=-1/2-yre+I yim;
Print[SetPrecision[rla,3],OutputForm[" "], ene[rla]//Re,
   OutputForm[" "],ene[rla]//Im,
   OutputForm[" "],d],
{rla,lamin,lamax,lastep}];
Clear[p,pol,t,s0,s,c,r,rt];
(* *)
(* Plotting ionization energy *)
(*SetDirectory["~/math/screened"];*)
redathe=Table[{la,Re[ene[la]]},{la,lamin,lamax,lastep}];
imdathe=Table[{la,Abs[Im[ene[la]]]},{la,lamin,lamax,lastep}];
fm={Interpolation[redat][la],Interpolation[imdat][la],
   Interpolation[redathe][la],Interpolation[imdathe][la]};
    (* plotting *)

<<Graphics`Legend`;
Clear[unit,leng,thic,xmin,xmax];
ticks[unit_,leng_,thic_][xmin_,xmax_]:=
   Block[{x,x0},
   x0=unit Ceiling[xmin/unit];
   Table[{x,x,{leng,0},{GrayLevel[0.],Thickness[thic]}},{x,x0,xmax,unit}]];
plt = 
 Plot[Release[fm],{la,lamin,lamax},
   Graphics`Legend`PlotLegend->{
StyleForm["\[Delta] = 0.88 (Re)","Section"],
StyleForm["\[Delta] = 0.88 (Im)","Section"],
StyleForm["He-like (Re)","Section"],
StyleForm["He-like (Im)","Section"]},
  Graphics`Legend`LegendPosition->
			{0.1,0.1},   
    TextStyle -> {FontFamily -> "Times", FontSize -> 24}, 
    PlotRange->{{lamin,lamax},{-0.1,0.5}},
    AxesLabel -> {\[Lambda], Subscript[E, I]}, 
    AxesStyle -> {Thickness[0.004]}, 
    AspectRatio->10/7.5,
    PlotStyle->{{Thickness[0.004],GrayLevel[0.5],RGBColor[1,0,0]},
{Thickness[0.004],GrayLevel[0.5],Dashing[{0.02,0.02}],RGBColor[1,0,0]},
{Thickness[0.004],RGBColor[0,0,1]},
{Thickness[0.004],Dashing[{0.02,0.02}],RGBColor[0,0,1]}},
    Ticks->{ticks[0.2,0.015,0.004],ticks[0.1,0.015,0.004]}
];
Display["plot3.gif",plt,"GIF",ImageSize->72 {15/2,10}];