(* Helium isoelectronic series *)
(* Finding resonances by variational method *)
(* With reading optimized a, b, c for low N *)
(* With calling FORTRAN subroutine to find eigenvalues *)

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

zlambda=1; (* 1 or 2, functions of Z or Lambda=1/Z *)
nm=25;
ndim=3; (* dimensionality - 3 or 5 or 7 ...*)
choice=1; (* 0 for i+j+k<=N, 1 for i+j^2+k^2<=N *)
ndig=96;
{tmin,tmax,tstep}={0.5,1,0.1};
inpf="opt5.dat";
ndig1=10; (* Number of printed digits *)
colmax0=2; (* expected optimum-convergence Pade column *)
colmax=1; (* max. Pade columns for printing *)
colmax1=3; (* max. Pade columns for plotting *)

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

nn=5; (* Degree of alg. appr. *)
ipts=22; (* Number of intermediate points for alg. interpolation *)
x3="   ";
SetOptions[$Output,PageWidth->Infinity];
outf="fix"<>ToString[nm]<>".dat";
eigdir="c:\\sergeev\\purdue\\fortran\\eigenval";
Clear[x];pri[x_]:=NumberForm[FortranForm[x//N],20];
{$MinPrecision,$MaxPrecision}={0,5000}; (* To avoid errors at running the following line *)
$MinPrecision=$MaxPrecision=ndig;
If[FileType[inpf]=!=File,  (* If the data file does not exist *)
 Print["(a,b,c)-data file does not exist"];
 Goto[99]
  ];
zcr=1/1.097660833;
(*zcrn=0.89158;*) (* for N = 3 *)
zcrn=0.89960; (* for N = 5 *)
{zcr,zcrn}=SetPrecision[{zcr,zcrn},ndig];
Clear[x];
chop[x_]=If[x==0,0,x];
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["Z-range: ",{zdt//Last,zdt//First}];
{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]]} };
Print["actual/data: ",MatrixForm[{{"N","Choice","Ndig"},{nm,choice,ndig},{nmdata,choicedata,ndigdata}}]];

os=OpenWrite[outf,PageWidth->Infinity,FormatType->OutputForm];
Write[os,nm,x3,choice,x3,ndig,x3,zlambda,x3,inpf];
Close[os];

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

Clear[ind,mmax];
mmm=hmapping[nm,choice,ind,mmax];
Print[Definition[mmm]];
{tmin,tmax,tstep}=Rationalize[{tmin,tmax,tstep}];
Do[ (* t = tmin - tmax, Z=t or Z=1/t *)
{tname,z,enscale}=If[zlambda==1,{"Z",t,1},{"1/Z",1/t,t^2}];
tpr=printF[t,8,4];
label=tname<>" = "<>tpr;
CellPrint[Cell[label,"Section",PageBreakAbove->True]];
z1=If[z<zcr&&z>0.85,z+zcrn-zcr,z]; (* shifting singularity to correct position *)
eig=alginterpolation[eguess,nn,z,intermediatePoints->ipts];
{a,b,c}={alginterpolation[apar,nn,z1,intermediatePoints->ipts],
         alginterpolation[bpar,nn,z1,intermediatePoints->ipts],
         alginterpolation[cpar,nn,z1,intermediatePoints->ipts]};
{apr,bpr,cpr}={printF[a,8,3],printF[b,8,3],printF[c,6,3]};
Print["{a,b,c} = ",{apr,bpr,cpr}];
time=Timing[
Clear[i1,j1,k1,i2,j2,k2];
{Tk,V0,V1,S}=hmatrix[a,b,c,i1,j1,k1,i2,j2,k2,ndim];
hintegrals[a,b,c,ndim,nm];
			]//First;
If[t==tmin,Print["time-integrals = ",time]];
A=B=T=Table[Null,{mmm},{mmm}];
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
B[[m1,m2]]=S;
,{m1,m2,mmm}],{m2,mmm}];
(* Orthogonalization *)
wvec=Table[Null,{n,mmm}];
Do[ (* Cycle through indexes of orthonormalized vectors *)
   T[[nn,nn]]=1;
   Do[wvec[[n]]=Sum[B[[nn,n1]] T[[n,n1]],{n1,n}],{n,nn-1}];					
   Do[T[[nn,n]]=-Sum[wvec[[n1]] T[[n1,n]],{n1,n,nn-1}],{n,nn-1}];					
   Do[wvec[[n]]=Sum[B[[n,n1]] T[[nn,n1]],{n1,n}]+
	  Sum[B[[n1,n]] T[[nn,n1]],{n1,n+1,nn}],{n,nn}];
   w=Sum[T[[nn,n]]wvec[[n]],{n,nn}];					
   norm=1/Sqrt[w];
   Do[T[[nn,n]]=norm T[[nn,n]],{n,nn}]
				,{nn,mmm}];
Clear[B];
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
A[[m1,m2]]=Tk-z V0+V1;
,{m1,m2,mmm}],{m2,mmm}];
(* Calculation of the matrix M *)
SetDirectory[eigdir];
matrs=OpenWrite["matr.dat",PageWidth->Infinity,FormatType->OutputForm];
Write[matrs,pri[z],x3,pri[eig//Re]," ",pri[eig//Im]];
Write[matrs,nm,x3,nmdata];
Do[Write[matrs,mmax[n]],{n,0,nm}];
time=Timing[
Do[
   Do[wvec[[k]]=
      Sum[T[[i,j]] A[[k,j]],{j,k}]+Sum[T[[i,j]] A[[j,k]],{j,k+1,i}],{k,i}];
   Do[result=Sum[wvec[[j]] T[[k,j]],{j,k}];
      Write[matrs,pri[result//Re]," ",pri[result//Im]];
      ,{k,i}],{i,mmm}];
			]//First;
Clear[A,T];
Close[matrs];
Run["eigenvc.bat"];
data=ReadList["eig.dat",Number,RecordSeparators->"\n",RecordLists->True];
ResetDirectory[];
If[t==tmin,Print["time-orthogonalization = ",time]];

(* Selecting the nearest eigenvalue *)
mseq=-1;
Do[
val0=enscale*(data[[n+2,1]]+I*data[[n+2,2]]);
enpr=printF[val0,16,12];
mseq++;
seq[mseq]=val0
     ,{n,0,nm}];
(* --- *)
Clear[x];
ser=seq[0]+Sum[(seq[n]-seq[n-1]) x^n,{n,1,mseq}];
tab=Table[" ",{l,mseq+2},{m,colmax+2}];
colm=Min[colmax0,mseq];
exact=Pade[ser,{x,0,mseq-colmax0,colm}]/.x->1;
enpr=printF[exact+0. I,16,12];
colmaxm=Max[colmax,colmax1];
Do[colm=Min[colmaxm,mseq-l];
	    Do[result=Pade[ser,{x,0,l,m}]/.x->1;
			pad[l,m]=result;
	   accpad[l,m]=-Log[10,Abs[result-exact]];
	If[l<=mseq&&m<=colmax,
		tab[[l+2,m+2]]=
        NumberForm[pad[l,m],{Infinity,ndig1},DigitBlock->3,
          NumberSeparator->" ",NumberSigns->{"-"," "}] ],{m,0,colm}];
If[l<=mseq,tab[[l+2,1]]=l];
	If[l<=colmax,tab[[1,l+2]]=l]
		,{l,0,mseq}];
Print[tab//MatrixForm];
accpr=printF[-Log[10,Abs[seq[mseq]-seq[mseq-1]]],6,2];
Do[lmax=If[m==colmax0,mseq-m-1(*to avoid infinity*),mseq-m];
   acclist[m]=Table[{l+m,accpad[l,m]},{l,0,lmax}],{m,0,colmax1}];
{$MinPrecision,$MaxPrecision}={0,5000};
MultipleListPlot@@Flatten[{Release@Table[acclist[m],{m,0,Min[colmax1,nm-1]}],
   PlotJoined->True,GridLines->Automatic},1];
$MinPrecision=$MaxPrecision=ndig;
os=OpenAppend[outf,PageWidth->Infinity,FormatType->OutputForm];
Print[tpr,x3,enpr,x3,accpr];
Write[os,tpr,x3,enpr,x3,accpr];
Unprotect[In,Out];Clear[In,Out];
Close[os]
	,{t,tmin,tmax,tstep}];
Label[99];
If[$MachineType!="PC",Exit[]];