Octonacci[k_]:=Module[{octo,i},
octo=Table[0,{i,1,k}];
octo[[1]]=1;
octo[[2]]=2;
Do[ octo[[i]]=octo[[i-2]]+2*octo[[i-1]],{i,3,k}];
octo];
locto=Octonacci[12];

locto

ClearAll[Mat,RawHamBloch,k,test,liens,cparmodcom]

liens=Import["liensPBC4v2.dat"];
cparmodcom=Import["PtsPBC4v2.dat"];
l=Length[liens];
ll=Length[cparmodcom];

cut=0.5;


RawHamBloch=Reap[Do[{rat=N[locto[[5]]*Sqrt[2]];dx=cparmodcom[[liens[[m,1]],1]]-cparmodcom[[liens[[m,2]],1]];dy=cparmodcom[[liens[[m,1]],2]]-cparmodcom[[liens[[m,2]],2]];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dx]>cut)&&(Abs[dy]>cut),(*Print["Cond 1 "," ",dx," ",dy," ",liens[[m,1]]," ",liens[[m,2]]];*)
Sow[{{liens[[m,1]],liens[[m,2]]},- Exp[I*rat*(kx*(dx-Sign[dx])+ky*(dy-Sign[dy]))]}];Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-I*rat*(kx*(dx-Sign[dx])+ky*(dy-Sign[dy]))]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dx]>cut)&&(Abs[dy]<cut),(*Print["Cond 2 "," ",dx," ",dy," ",liens[[m,1]]," ",liens[[m,2]]];*)
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[I*rat*(kx*(dx-Sign[dx])+ky*(dy))]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-I*rat*(kx*(dx-Sign[dx])+ky*(dy))]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dy]>cut)&&(Abs[dx]<cut),(*Print["Cond 3 "," ",dx," ",dy," ",liens[[m,1]]," ",liens[[m,2]]];*)
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[I*rat*(kx*(dx)+ky*(dy-Sign[dy]))]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-I*rat*(kx*(dx)+ky*(dy-Sign[dy]))]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dy]<cut)&&(Abs[dx]<cut),(*Print["Cond 4 "," ",dx," ",dy," ",liens[[m,1]]," ",liens[[m,2]]];*)
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[I*rat*(kx*(dx)+ky*(dy))]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-I*rat*(kx*(dx)+ky*(dy))]}];]
},{m,1,l}]][[2,1]];

tmp=Split[Sort[RawHamBloch],#1[[1]]==#2[[1]]&];
RawHamBloch=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];

(***** Test convergence des résultats ******)
rulef={kx->0.,ky->0.};
ham=SparseArray[RawHamBloch/.rulef];

rulef2={kx->N[Pi/(3*70*Sqrt[2])],ky->0.};
ham2=SparseArray[RawHamBloch/.rulef];

Sort[Eigenvalues[ham,3,Method->{"Arnoldi",MaxIterations->1000}]]

Sort[Eigenvalues[ham,3,Method->{"Arnoldi",MaxIterations->4000}]]

{-4.22 1690550775537`,4.219022342720782`,4.21914987375085`}

Sort[Eigenvalues[ham]]

Sort[Eigenvalues[ham,10,Method->{"Arnoldi",MaxIterations->1000}]]

{-4.22 16905507749924`,-4.216442578873379`,-4.215974714907718`,-4.215739378039369`,-4.215031267559585`,-4.211209391398773`,4.218160598903179`,4.218592268291222`,4.219022353455503`,4.219149873953131`}

Sort[Eigenvalues[ham,10,Method->{"Arnoldi",MaxIterations->4000}]]

{-4.221 690550775144`,-4.216442578882125`,-4.2159747180404805`,-4.215739379510237`,-4.2150312682140525`,-4.211213852955974`,4.218160598904114`,4.218592268292102`,4.219022353456116`,4.219149873953587`}

Sort[Eigenvalues[ham2,10,Method->{"Arnoldi",MaxIterations->1000}]]

{-4.22 16905507749924`,-4.216442578873379`,-4.215974714907718`,-4.215739378039369`,-4.215031267559585`,-4.211209391398773`,4.218160598903179`,4.218592268291222`,4.219022353455503`,4.219149873953131`}

(***** Choix fenêtre de fit *****)

locto[[5]]

N[Pi/(3*29*Sqrt[2])]

sp5xv1=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*N[Pi/(6*70*Sqrt[2])],ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(6*70*Sqrt[2])],Sort[Eigenvalues[ham,10]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],n];

sp5xv2=Monitor[Reap[
For[n=1,n≤9,n++,
rulef={kx->(n-5)*N[Pi/(16*70*Sqrt[2])],ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-5)*N[Pi/(16*70*Sqrt[2])],Sort[Eigenvalues[ham,10]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],n];

sp5x=Join[sp5xv2,sp5xv1];

Show[Plot[-4.221690548524171`+1.443000089518171` (x^2)+3.597617557352529`*^-7 (x^3),{x,-0.006,0.006}],ListPlot[Take[Transpose[{Transpose[sp5xu][[1]],Transpose[sp5xu][[2]]}],{3,11}]]]

Fit[Transpose[{Transpose[sp5xu][[1]],Transpose[sp5xu][[2]]}],{1,x^2,x^3},x]

Fit[Take[Transpose[{Transpose[sp5xu][[1]],Transpose[sp5xu][[2]]}],{2,12}],{1,x^2,x^3},x]

Fit[Take[Transpose[{Transpose[sp5xu][[1]],Transpose[sp5xu][[2]]}],{3,11}],{1,x^2,x^3},x]

sp5xu=Union[sp5x];

sp5yv1=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={ky->(n-3)*N[Pi/(6*70*Sqrt[2])],kx->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(6*70*Sqrt[2])],Sort[Eigenvalues[ham,10]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],n];
sp5yv2=Monitor[Reap[
For[n=1,n≤9,n++,
rulef={ky->(n-5)*N[Pi/(16*70*Sqrt[2])],kx->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-5)*N[Pi/(16*70*Sqrt[2])],Sort[Eigenvalues[ham,10]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],n];
sp5y=Union[Join[sp5yv1,sp5yv2]];

ListPlot[Take[Transpose[{Transpose[sp5y][[1]],Transpose[sp5y][[2]]}],{3,11}]]

Fit[Take[Transpose[{Transpose[sp5y][[1]],Transpose[sp5y][[2]]}],{1,13}],{1,x^2,x^3},x]

Fit[Take[Transpose[{Transpose[sp5y][[1]],Transpose[sp5y][[2]]}],{3,11}],{1,x^2,x^3},x]

Fit[Take[Transpose[{Transpose[sp5y][[1]],Transpose[sp5y][[2]]}],{4,10}],{1,x^2,x^3},x]

sp5xyv1=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={ky->(n-3)*N[Pi/(6*70*Sqrt[2])],kx->(n-3)*N[Pi/(6*70*Sqrt[2])]};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(6*70*Sqrt[2])],Sort[Eigenvalues[ham,10]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],n];
sp5xyv2=Monitor[Reap[
For[n=1,n≤9,n++,
rulef={ky->(n-5)*N[Pi/(16*70*Sqrt[2])],kx->(n-5)*N[Pi/(16*70*Sqrt[2])]};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-5)*N[Pi/(16*70*Sqrt[2])],Sort[Eigenvalues[ham,10]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],n];
sp5xy=Union[Join[sp5xyv1,sp5xyv2]];

ListPlot[Take[Transpose[{Transpose[sp5xy][[1]],Transpose[sp5xy][[2]]}],{3,11}]]

+1.414363189469766`

Fit[Take[Transpose[{Transpose[sp5xy][[1]],Transpose[sp5xy][[2]]}],{3,11}],{1,x^2,x^3},x]

Fit[Take[Transpose[{Transpose[sp5xy][[1]],Transpose[sp5xy][[2]]}],{4,10}],{1,x^2,x^3},x]

alpha5x=2*1.443000089518171`
alpha5y=2*1.4431009749827806`
alpha5xy=2.8859127047383755`-1.443000089518171`-1.4431009749827806`

2.886000179036342`

2.886201949965561`

-0.0001 8835976257602915`

minv5=Sqrt[alpha5x*alpha5y-alpha5xy*alpha5xy]

(**** Approximant k=4 *****)

sp4xv1=Monitor[Reap[
For[n=1,n≤21,n++,
rulef={kx->(n-11)*N[Pi/(4*10*29*Sqrt[2])],ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-11)*N[Pi/(4*10*29*Sqrt[2])],Sort[Eigenvalues[ham,10]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],n];


sp4yv1=Monitor[Reap[
For[n=1,n≤21,n++,
rulef={ky->(n-11)*N[Pi/(40*29*Sqrt[2])],kx->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-11)*N[Pi/(40*29*Sqrt[2])],Sort[Eigenvalues[ham,10]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],n];

sp4xyv1=Monitor[Reap[
For[n=1,n≤21,n++,
rulef={ky->(n-11)*N[Pi/(40*29*Sqrt[2])],kx->(n-11)*N[Pi/(40*29*Sqrt[2])]};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-11)*N[Pi/(40*29*Sqrt[2])],Sort[Eigenvalues[ham,10]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],n];

ListPlot[Take[Transpose[{Transpose[sp4xv1][[1]],Transpose[sp4xv1][[2]]}],{7,15}]]
Fit[Take[Transpose[{Transpose[sp4xv1][[1]],Transpose[sp4xv1][[2]]}],{6,15}],{1,x^2,x^3},x]

Fit[Take[Transpose[{Transpose[sp4xv1][[1]],Transpose[sp4xv1][[2]]}],{7,15}],{1,x^2,x^3},x]

Fit[Take[Transpose[{Transpose[sp4yv1][[1]],Transpose[sp4yv1][[2]]}],{7,15}],{1,x^2,x^3},x]

Fit[Take[Transpose[{Transpose[sp4xyv1][[1]],Transpose[sp4xyv1][[2]]}],{7,15}],{1,x^2,x^3},x]

alpha4x=2*1.4586853102365855`;
alpha4y=2*1.4586853088986131`;
alpha4xy=2.917750400003968`-1.4586853102365855`-1.4586853088986131`

minv4=Sqrt[alpha4x*alpha4y-alpha4xy*alpha4xy]