(* Calculating of the energy (including resonances) for helium isoelectronic 
    series by finding a stationary point of the variational functional *)
(* With advanced differentiation over a, b, c *)
(* With optimized performance *)
(* Using algebraic extrapolation for guessing *)
(* As a function of Z, without re-scaling *)
(* Finding Zcr where E+Z^2/2 = 0 or Z* where Im E = 0 *)

(********** INPUT STARTS **********)

ibreak=2; (* 1 for Re E+Z^2/2 > 0, 2 for |Im E| > 0 *)
{nmin,nmax}={0,0};
z0=0.89501;
zstep=0.03; (* The critical Z should be in the interval (z0-zstep,z0+zstep) *)
ipts=9; (* Number of intermediate points for alg. interpolation *)
maxdegree=3; (* Max. degree of algebraic interpolation *)
choice=0; (* 0 for i+j+k<=N, 1 for i+j^2+k^2<=N *)
ndig=96; (* Precision of calculations *)

(**********  INPUT ENDS  **********)

epsim=10.^-10; (* Allowance of im. parts *)
ndig1=Min[ndig-10,22]; (* AccuracyGoal in FindRoot *)
ndig2=10; (* Accuracy of finding the cr. point *)
ndim=3; (* dimensionality - 3 or 5 or 7 ...*)
miter=100; (* max n. of iterations *)
x=" ";
x3="   ";
outf1=Switch[ibreak,1,"z-cr.dat",2,"z-ast.dat"];
{$MinPrecision,$MaxPrecision}={0,5000}; (* To avoid errors at running the following line *)
$MinPrecision=$MaxPrecision=ndig;
Off[Precision::"precsm"];
Off[Solve::"precsm"];
Clear[t];
chop[t_]=If[t==0,0,t];
{z0,zstep}=Rationalize[{z0,zstep}];
Print["{Nrange,choice,precision} = ",{{nmin,nmax},If[choice==0," i+j+k<=N "," i+j^2+k^2<=N "],ndig}];

p=ToFileName["..","packages"];
If[!MemberQ[$Path,p],$Path=Append[$Path,p]];
<<hylleraa`;
<<alginter`;
<<printing` (* no ";" *)

Clear[a,b,c,i,j,k,i1,j1,k1,i2,j2,k2,int1,int2,ints,int1a,int2a,intsa,int1b,int2b,intsb,int1c,int2c,intsc];
{T,V0,V1,S}=hmatrix[a,b,c,i1,j1,k1,i2,j2,k2,ndim];
(* Differentiating the matrix *)
{T,V0,V1,S}={T,V0,V1,S}/.{
   hylleraas`private`int1[i_,j_,k_]->int1[i,j,k][a,b,c],
   hylleraas`private`int2[i_,j_,k_]->int2[i,j,k][a,b,c],
   hylleraas`private`ints[i_,j_,k_]->ints[i,j,k][a,b,c]
		};
{Ta,V0a,V1a,Sa}=D[{T,V0,V1,S},a]/.{
		Derivative[1,0,0][int1[i_,j_,k_]][a,b,c]->int1a[i,j,k],
		Derivative[1,0,0][int2[i_,j_,k_]][a,b,c]->int2a[i,j,k],
		Derivative[1,0,0][ints[i_,j_,k_]][a,b,c]->intsa[i,j,k],
		int1[i_,j_,k_][a,b,c]->int1[i,j,k],
		int2[i_,j_,k_][a,b,c]->int2[i,j,k],
		ints[i_,j_,k_][a,b,c]->ints[i,j,k]
		};
{Tb,V0b,V1b,Sb}=D[{T,V0,V1,S},b]/.{
		Derivative[0,1,0][int1[i_,j_,k_]][a,b,c]->int1b[i,j,k],
		Derivative[0,1,0][int2[i_,j_,k_]][a,b,c]->int2b[i,j,k],
		Derivative[0,1,0][ints[i_,j_,k_]][a,b,c]->intsb[i,j,k],
		int1[i_,j_,k_][a,b,c]->int1[i,j,k],
		int2[i_,j_,k_][a,b,c]->int2[i,j,k],
		ints[i_,j_,k_][a,b,c]->ints[i,j,k]
		};
{Tc,V0c,V1c,Sc}=D[{T,V0,V1,S},c]/.{
		Derivative[0,0,1][int1[i_,j_,k_]][a,b,c]->int1c[i,j,k],
		Derivative[0,0,1][int2[i_,j_,k_]][a,b,c]->int2c[i,j,k],
		Derivative[0,0,1][ints[i_,j_,k_]][a,b,c]->intsc[i,j,k],
		int1[i_,j_,k_][a,b,c]->int1[i,j,k],
		int2[i_,j_,k_][a,b,c]->int2[i,j,k],
		ints[i_,j_,k_][a,b,c]->ints[i,j,k]
		};
{T,V0,V1,S}={T,V0,V1,S}/.{
			int1[i_,j_,k_][a,b,c]->int1[i,j,k],
			int2[i_,j_,k_][a,b,c]->int2[i,j,k],
			ints[i_,j_,k_][a,b,c]->ints[i,j,k]
		};
Do[ (* nm = nmin, ..., nmax *)
Print[nm//Definition];
outf="t"<>ToString[nm]<>".dat";
If[FileType[outf]=!=File, (* If the data file does not exist *)
 Print["Input data file does not exist!";Continue[]]];
inpf=outf;
data=ReadList[inpf,Number,RecordSeparators->"\n",RecordLists->True];
data=ReadList[inpf,Number,RecordSeparators->"\n",RecordLists->True];
(* One line does not work correctly for unknown reason.
   It is related somehow to setting $MinPrecision and $MaxPrecision. *)
{nmdata,choicedata,ndigdata}=data//First;
data=Drop[data,1];
dt=Transpose[data];
zdt=dt[[1]]//Rationalize; (* Rounding Z to exact number *)
dt=SetPrecision[dt,ndig];
Print["Using the previous output file ",outf];
Print["Z-previous: ",printF[zdt//First,8,4],"  - ",printF[zdt//Last,8,4],
    "   (",Length[zdt]," data points)"];
{apar,bpar,cpar,eguess}=Transpose/@{{zdt,dt[[2]]+I chop/@dt[[3]]},
    {zdt,dt[[4]]+I chop/@dt[[5]]},
    {zdt,dt[[6]]+I chop/@dt[[7]]},
    {zdt,dt[[8]]+I chop/@dt[[9]]} };
(**)  
Clear[ind,mmax];
mm=hmapping[nm,choice,ind,mmax];
Print[Definition[mm]];
L=M=Table[Null,{mm},{mm}];
(* *)
Clear[a,b,c];
time0=Timing[
hintegrals[a,b,c,ndim,nm];
(* Diffirentiating integrals *)
nmax=2 nm;
nmadd=2 (ndim-3);
nmax3=nmax+nmadd+3;
Do[
  Do[namax=n-nc;namax1=Floor[namax/2];
    Do[nb=namax-na;
       integ1[na-1,nb-1,nc-1]=integ2[nb-1,na-1,nc-1]=hylleraas`private`int1[na-1,nb-1,nc-1];
       integ2[na-1,nb-1,nc-1]=integ1[nb-1,na-1,nc-1]=hylleraas`private`int2[na-1,nb-1,nc-1];
       integs[na-1,nb-1,nc-1]=integs[nb-1,na-1,nc-1]=hylleraas`private`ints[na-1,nb-1,nc-1];
       integ1a[na-1,nb-1,nc-1]=integ2a[nb-1,na-1,nc-1]=
              D[hylleraas`private`int1[na-1,nb-1,nc-1],a]//Together;
       integ2a[na-1,nb-1,nc-1]=integ1a[nb-1,na-1,nc-1]=
              D[hylleraas`private`int2[na-1,nb-1,nc-1],a]//Together;
       integsa[na-1,nb-1,nc-1]=integsa[nb-1,na-1,nc-1]=
              D[hylleraas`private`ints[na-1,nb-1,nc-1],a]//Together;
       integ1b[na-1,nb-1,nc-1]=integ2b[nb-1,na-1,nc-1]=
              D[hylleraas`private`int1[na-1,nb-1,nc-1],b]//Together;
       integ2b[na-1,nb-1,nc-1]=integ1b[nb-1,na-1,nc-1]=
              D[hylleraas`private`int2[na-1,nb-1,nc-1],b]//Together;
       integsb[na-1,nb-1,nc-1]=integsb[nb-1,na-1,nc-1]=
              D[hylleraas`private`ints[na-1,nb-1,nc-1],b]//Together;
       integ1c[na-1,nb-1,nc-1]=integ2c[nb-1,na-1,nc-1]=
              D[hylleraas`private`int1[na-1,nb-1,nc-1],c]//Together;
       integ2c[na-1,nb-1,nc-1]=integ1c[nb-1,na-1,nc-1]=
              D[hylleraas`private`int2[na-1,nb-1,nc-1],c]//Together;
       integsc[na-1,nb-1,nc-1]=integsc[nb-1,na-1,nc-1]=
              D[hylleraas`private`ints[na-1,nb-1,nc-1],c]//Together;
      ,{na,0,namax1}],{nc,0,n}],{n,0,nmax3}];
]//First;Print["time-matrix-general form = ",time0];
(* *)
z=z0;
zst=zstep;
While[zst>10.^(-ndig2), (* Z-loop starts here *)
(* Guessing *)
mdata=Length[eguess];
Clear[n];
nn=Max[Floor[n/.Solve[(n+1)(n+2)/2-1==mdata,n]]];
nn=Min[nn,maxdegree]; (* up to max. degree *)
Print["Z = ",printF[z,ndig2+4,ndig2],"   Zstep = ",printF[zst,ndig2+4,ndig2],"     degree of alg. f. = ",nn];
eig=alginterpolation[eguess,nn,z,intermediatePoints->ipts];
abc1=  {alginterpolation[apar,nn,z,intermediatePoints->ipts],
       alginterpolation[bpar,nn,z,intermediatePoints->ipts],
       alginterpolation[cpar,nn,z,intermediatePoints->ipts]};
abc2= {alginterpolation[apar,nn-1,z,intermediatePoints->ipts],
       alginterpolation[bpar,nn-1,z,intermediatePoints->ipts],
       alginterpolation[cpar,nn-1,z,intermediatePoints->ipts]}
      +10^-11 (1+I) {1,1,1};
guess={{atrial,{abc1[[1]],abc2[[1]]}},
       {btrial,{abc1[[2]],abc2[[2]]}},
       {ctrial,{abc1[[3]],abc2[[3]]}}};
(******************************* START find-root ******************************)
dampf=1;
iter=0;
abciter={};
rt=FindRoot[
iter++;
{a,b,c}=
   SetPrecision[{atrial,btrial,ctrial},ndig]//N; (* increase of precision *)
(* finding the eigenvalue *)
time1=Timing[
Clear[int1,int2,ints,int1a,int2a,intsa,int1b,int2b,intsb,int1c,int2c,intsc,na,nb,nc];
int1[na_,nb_,nc_]:=int1[na,nb,nc]=int2[nb,na,nc]=integ1[na,nb,nc];
int2[na_,nb_,nc_]:=int2[na,nb,nc]=int1[nb,na,nc]=integ2[na,nb,nc];
ints[na_,nb_,nc_]:=ints[na,nb,nc]=ints[nb,na,nc]=integs[na,nb,nc];
int1a[na_,nb_,nc_]:=int1a[na,nb,nc]=int2a[nb,na,nc]=integ1a[na,nb,nc];
int2a[na_,nb_,nc_]:=int2a[na,nb,nc]=int1a[nb,na,nc]=integ2a[na,nb,nc];
intsa[na_,nb_,nc_]:=intsa[na,nb,nc]=intsa[nb,na,nc]=integsa[na,nb,nc];
int1b[na_,nb_,nc_]:=int1b[na,nb,nc]=int2b[nb,na,nc]=integ1b[na,nb,nc];
int2b[na_,nb_,nc_]:=int2b[na,nb,nc]=int1b[nb,na,nc]=integ2b[na,nb,nc];
intsb[na_,nb_,nc_]:=intsb[na,nb,nc]=intsb[nb,na,nc]=integsb[na,nb,nc];
int1c[na_,nb_,nc_]:=int1c[na,nb,nc]=int2c[nb,na,nc]=integ1c[na,nb,nc];
int2c[na_,nb_,nc_]:=int2c[na,nb,nc]=int1c[nb,na,nc]=integ2c[na,nb,nc];
intsc[na_,nb_,nc_]:=intsc[na,nb,nc]=intsc[nb,na,nc]=integsc[na,nb,nc];
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
L[[m1,m2]]=L[[m2,m1]]=T-z V0+V1;
M[[m1,m2]]=M[[m2,m1]]=S;
,{m1,m2,mm}],{m2,mm}];
A=Inverse[M].L;
			]//First;(*Print["time-matrix = ",time1];*)
time2=Timing[
{val0,vec0}=
              Sort[Transpose[Eigensystem[A]],
                  Abs[#1[[1]]-eig]<Abs[#2[[1]]-eig]&][[1]];
			]//First;(*Print["time-eigenvalues = ",time2];*)
(* finding derivative of the eigenvalue d/da *)
time3=Timing[
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
L[[m1,m2]]=L[[m2,m1]]=Ta-z V0a+V1a;
M[[m1,m2]]=M[[m2,m1]]=Sa;
,{m1,m2,mm}],{m2,mm}];
vala=(vec0.L.vec0)/(vec0.M.vec0)-val0;
(* finding derivative of the eigenvalue d/db *)
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
L[[m1,m2]]=L[[m2,m1]]=Tb-z V0b+V1b;
M[[m1,m2]]=M[[m2,m1]]=Sb;
,{m1,m2,mm}],{m2,mm}];
valb=(vec0.L.vec0)/(vec0.M.vec0)-val0;
(* finding derivative of the eigenvalue d/dc *)
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
L[[m1,m2]]=L[[m2,m1]]=Tc-z V0c+V1c;
M[[m1,m2]]=M[[m2,m1]]=Sc;
,{m1,m2,mm}],{m2,mm}];
valc=(vec0.L.vec0)/(vec0.M.vec0)-val0;
			]//First;
delt={vala,valb,valc};
deltabs=Plus@@Abs/@delt;
deltpr=PaddedForm[deltabs//FortranForm,{3,1},ExponentFunction->Identity];
abciter=Append[abciter,{deltabs,{atrial,btrial,ctrial},val0}];
If[iter==1&&z==zstart,Print["Time (one step): matrix, eigenv., derivatives:  ",time1,", ",time2,", ",time3]];
Print[iter,"  --- ",deltpr,"  ",{printF[atrial,9,6],printF[btrial,9,6],printF[ctrial,9,6]}
     ];delt=={0,0,0},
      Sequence@@guess//Release,
WorkingPrecision->ndig,AccuracyGoal->ndig1,MaxIterations->miter,DampingFactor->dampf];
(******************************* END find-root ******************************)
{delt,abc0,eig0}=Sort[abciter,#1[[1]]<#2[[1]]&][[1]];
{a0,b0,c0}=abc0;
zpr=printF[z,ndig2+2,ndig2-2];
{apr,bpr,cpr}={printF[a0+0. I,ndig1+2,ndig1-2],
 printF[b0+0. I,ndig1+2,ndig1-2],
 printF[c0+0. I,ndig1+2,ndig1-2]};
enpr=printF[eig0+0. I,ndig1+4,ndig1];
deltpr=PaddedForm[delt//FortranForm,{3,1},ExponentFunction->Identity];
Print["===========   Z = ",zpr];
Print["Non-linear parameters:  ",{apr,bpr,cpr}];
Print["Energy:  ",enpr];
Print["L.h.s. of the equations (should be zero), n. of iter.:  ",deltpr,", ",iter];
If[iter>=miter,
 Print["Number of iterations exceeded the limit of ",miter," iterations."];
 Break[] ];
os=OpenAppend[outf,PageWidth->Infinity,FormatType->OutputForm];
Write[os,zpr,x3,apr,x,bpr,x,cpr,x3,enpr,x3,deltpr];
Close[os];
{apar,bpar,cpar,eguess}=
  {Append[apar,{z,a0}],Append[bpar,{z,b0}],Append[cpar,{z,c0}],
   Append[eguess,{z,eig0}]};
break=Switch[ibreak,
      1,Re[eig0]+z^2/2>0,
      2,Max[Abs[Im[a0]],Abs[Im[b0]],Abs[Im[c0]],Abs[Im[eig0]]]>epsim];
z=If[break,z+zst,z-zst];
zst=zst/2
]; (* End of Z While-loop *)
os=OpenAppend[outf1,PageWidth->Infinity,FormatType->OutputForm];
Write[os,nm,x3,zpr,x3,apr,x,bpr,x,cpr,x3,enpr,x3,deltpr];
Close[os];
,{nm,nmin,nmax}];
If[$MachineType!="PC",Exit[]];