(* Interpolation of energy levels and adding additional points
   for subsequent plotting *)

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

nm=12; (* N *)
mlevels=10; (* Max. number of treated levels *)
ndeg=3; (* Degree of alg. interp. *)
ipts=9; (* N. of intermed. points for calc. of interpolation *)
zrange=Flatten[{
Range[0.6,0.8,0.005],Range[0.802,0.85,0.002],
Range[0.851,0.9,0.001],Range[0.9005,0.91,0.0005],Range[0.911,0.92,0.001],
Range[0.925,1,0.005]
       }];
ndig=16; (* Precision of calculations *)

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

zrange=zrange//Rationalize;
{$MinPrecision,$MaxPrecision}={0,5000}; (* To avoid errors at running the following line *)
$MinPrecision=$MaxPrecision=ndig;
inpf="fixabc"<>ToString[nm]<>".dat";
data=ReadList[inpf,Number,RecordSeparators->"\n",RecordLists->True];
{nmdata,choicedata,ndigdata,a,b,c}=data//First;
data=Drop[data,1];
dt=Transpose[data];
{zrange,dt}=SetPrecision[{zrange,dt},ndig];
zdt=dt[[1]]//Rationalize; (* Rounding Z to exact number *)
mlev=Min[mlevels,Length[data[[1]]]-1];
Clear[lev];
Release[lev/@Range[mlev]]=Transpose/@Table[{zdt,dt[[n+1]]},{n,mlev}];
outf="f0abc"<>ToString[nm]<>".dat";
Close[outf];
os=OpenWrite[outf,PageWidth->Infinity,FormatType->OutputForm];

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


mz=Length[zrange];
mzdt=Length[zdt];
Do[z=zrange[[nz]];
   nn=ndeg;
   Do[If[zdt[[nzdt]]==z,nn=0;Break[]],{nzdt,mzdt}];
   line=printF[z,9,5]<>"  ";
   Print["Z =",line];
   Do[line=line<>" "<>printF[
         Re[alginterpolation[lev[nlev],nn,z,intermediatePoints->ipts]]
                            ,12,8],{nlev,mlev}];
   Write[os,line]
   ,{nz,mz}];
Close[os];