(* Helium isoelectronic series *)
(* Finding resonances by variational method *)
(* With fixed real a,b,c for all Z (for testing) *)

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

zlambda=1; (* 1 or 2, functions of Z or Lambda=1/Z *)
nm=12;
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;
mlevels=100; (* Max. number of levels *)
{afix,bfix,cfix}={0.1,0.8,0};
{tmin,tmax,tstep}={0.6,1,0.005};

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

x3="   ";
SetOptions[$Output,PageWidth->Infinity];
outf="fixabc"<>ToString[nm]<>".dat";
{$MinPrecision,$MaxPrecision}={0,5000}; (* To avoid errors at running the following line *)
$MinPrecision=$MaxPrecision=ndig;
{afix,bfix,cfix}=SetPrecision[{afix,bfix,cfix},ndig];
Clear[x];
chop[x_]=If[x==0,0,x];
Print[MatrixForm[{{"N","Choice","Ndig"},{nm,choice,ndig}}]];

os=OpenWrite[outf,PageWidth->Infinity,FormatType->OutputForm];
Write[os,nm,x3,choice,x3,ndig,x3,afix,x3,bfix,x3,cfix];
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}];
A=B=M0=T=Table[Null,{mmm},{mmm}];
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;
Print[label];
{a,b,c}={afix,bfix,cfix};
{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]];
time=Timing[
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
A[[m1,m2]]=Tk-z V0+V1;
B[[m1,m2]]=S;
,{m1,m2,mmm}],{m2,mmm}];
			]//First;
If[t==tmin,Print["time-matrixes = ",time]];
(* Orthogonalization *)
time=Timing[
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}];
(* Calculation of the matrix M *)
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}];
      M0[[i,k]]=result//N;
      ,{k,i}],{i,mmm}];
			]//First;
If[t==tmin,Print["time-orthogonalization = ",time]];
(* Selecting the nearest eigenvalue *)
mm=mmax[nm];
M=Table[Null,{mm},{mm}];
Do[Do[ M[[i,k]]=M[[k,i]]=M0[[i,k]],{k,i}],{i,mm}];
{$MinPrecision,$MaxPrecision}={0,5000};
eigs=Eigenvalues[M//N];
$MinPrecision=$MaxPrecision=ndig;
eigsort=Sort[eigs];
Clear[enpr];
nlevels=Min[mlevels,mm];
pri="";
Do[pri=pri<>" "<>printF[enscale*(-eigsort[[i]]-z^2/2),16,12],{i,nlevels}]; (* Printing ionization energy *)
(* --- *)
os=OpenAppend[outf,PageWidth->Infinity,FormatType->OutputForm];
Write[os,tpr,x3,pri];
Close[os]
	,{t,tmin,tmax,tstep}];
Label[99];
If[$MachineType!="PC",Exit[]];