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]

(*********** Matrice définie avec SparseArray ****************)
liens=Import["liensPBC3.dat"];
cparmodcom=Import["PtsPBC3.dat"];
l=Length[liens];
ll=Length[cparmodcom];

cut=0.5;


RawHamBloch=Reap[Do[{rat=N[locto[[7]]*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]}];

ll

ListPlot[Transpose[{Transpose[sp3xv3][[1]],Transpose[sp3xv3][[2]]}]]

ListPlot[Transpose[{Transpose[sp6xv3][[1]],Transpose[sp6xv3][[2]]}]]

N[Pi/(169*Sqrt[2])]

alphax6v5=2*1.4082910300528955`;
alphay6v5=2*1.4097308131227646`;
alphaxyv5=2.793556709771991`-1.4082910300528955`-1.4097308131227646`

minv6v5=Sqrt[alphax6v5*alphay6v5-alphaxyv5*alphaxyv5]

Pi*2.8179152742563778`

Fit[Transpose[{Transpose[sp6xv6][[1]],Transpose[sp6xv6][[2]]}],{1,x^2},x]

Clear[sp6xv7,sp6yv7,sp6xyv7]

sp6xv11=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)/200.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)/200.,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[Transpose[{Transpose[sp6xv12][[1]],Transpose[sp6xv12][[2]]}]]

sp6xv12=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)/20.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)/20.,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];

N[Pi/(169*Sqrt[2])]

sp6xv13=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)/80.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)/80.,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[Transpose[{Transpose[sp6xv13][[1]],Transpose[sp6xv13][[2]]}]]

sp6xv14=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*13/4000.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*13/4000.,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[Transpose[{Transpose[sp6xv15][[1]],Transpose[sp6xv15][[2]]}]]

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

sp6xv15=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*13/6000.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*13/6000.,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];

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

sp6xv16=Join[sp6xv15,sp6xv14];

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

sp6xv17=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*13/8000.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*13/8000.,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];

N[26/8000]

ListPlot[Transpose[{Transpose[sp6xv17][[1]],Transpose[sp6xv17][[2]]}]]

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

sp6xv18=Monitor[Reap[
For[n=1,n≤4,n++,
rulef={kx->n*13/32000.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{n*13/32000.,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[Transpose[{Transpose[sp6xv18][[1]],Transpose[sp6xv18][[2]]}]]

sp6fil=Take[sp6,{3,6}]
Fit[Transpose[{Transpose[sp6fil][[1]],Transpose[sp6fil][[2]]}],{1,x,x^2,x^3},x]

1.425*2

sp6x=Join[sp6xv14,sp6xv15,sp6xv17,sp6xv18];

Show[ListPlot[Transpose[{Transpose[spi][[1]],Transpose[spi][[2]]}]],ListPlot[Transpose[{Transpose[sp6xf][[1]],Transpose[sp6xf][[2]]}],PlotStyle->Red]]

Export["Sp6xBlochv2.dat",Union[N[sp6x]]];

spi2=Import["Sp6xBlochv2.dat"];

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

sp6xf=Take[Sort[Union[sp6x]],{3,12}];

Length[spi2]

N[13/4000]

0.003 25`

Union[Transpose[{Transpose[spi2][[1]],Transpose[spi2][[2]]}],Transpose[{-Transpose[spi2][[1]],Transpose[spi2][[2]]}]]

Length[spicompl]

ListPlot[Take[spicompl,{5,13}]]

Show[{ListPlot[Take[spicompl,{4,14}]],Plot[-4.221696286655369`+1.4268397249646436` (x^2)+0.009284676785805626` (x^3)-88.84770861592226` (x^4),{x,-0.006,0.006}]}]

Length[sp6c]

ListPlot[Transpose[{Transpose[sp6c][[1]],Transpose[sp6c][[2]]}]]

sp6cc

Length[Union[sp6x]]

sp6xv19=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*325/200000.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*325/200000.,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];

N[2*325/200000]

sp6yv1=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={ky->(n-3)*13/8000.,kx->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*13/8000.,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];

sp6yv2=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={ky->(n-3)*13/6000.,kx->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*13/6000.,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];
sp6comp2=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*13/6000.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*13/6000.,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];

sp6y=Union[Join[sp6yv3,sp6yv2,sp6yv1]];

sp6comp=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*13/8000.,ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*13/8000.,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];

sp6c=Union[sp6comp,sp6comp2,sp6comp3];

ListPlot[Transpose[{Transpose[sp6y][[1]],Transpose[sp6y][[2]]}]]

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

Length[sp6y]

Fit[Take[Transpose[{Transpose[sp6y][[1]],Transpose[sp6y][[2]]}],{3,9}],{1,x^2,x^3},x]

ListPlot[Take[Transpose[{Transpose[sp6y][[1]],Transpose[sp6y][[2]]}],{3,9}]]

sp6xyv1=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*13/6000.,ky->(n-3)*13/6000.};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*13/6000.,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];
sp6xyv2=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*13/16000.,ky->(n-3)*13/16000.};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*13/16000.,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];

sp6xy=Union[Join[sp6xyv2,sp6xyv1]];

Length[sp6xy]

ListPlot[Take[Transpose[{Transpose[sp6xy][[1]],Transpose[sp6xy][[2]]}],{2,8}]]

Fit[Take[Transpose[{Transpose[sp6xy][[1]],Transpose[sp6xy][[2]]}],{2,8}],{1,x^2,x^3},x]

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

alpha6x=2*1.4265963300109266`;
alpha6y=2*1.4263898626442608`;

alpha6xy=2.8526222672928676`-1.4265963300109266`-1.4263898626442608`;

minv6=Sqrt[alpha6x*alpha6y-alpha6xy*alpha6xy]

(****** Convergence des calculs numériques ******)

rulef={kx->0.,ky->0.};
ham=SparseArray[RawHamBloch/.rulef];

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

{-4.221617639351884`, 4.220958228550611`,4.221173568294972`}

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

{-4.221 6962866595`,-4.220748952466333`,-4.220668364527309`,-4.2197196167913384`,-4.21911851072796`,4.2188682476198505`,4.2194538271195725`,4.2207723295795025`,4.221144720585183`,4.221232655908404`}

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

{-4.2216962856581866`,4 .221147299948926`,4.221237243696406`}

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

{-4.221696286973762`,- 4.220821717346861`,-4.220705134505429`,-4.220665248731077`,-4.220583377192047`,-4.219887511750188`,4.221121522372185`,4.221142273622645`,4.221220981651449`,4.221288366447012`}

(* Conclusion : on choisit, pour k=6, de faire les calculs en prenant 10 vap, 1000 itérations*)

N[13/4000]

rulef2={kx->0.0032500000000000003`,ky->0.};
ham2=SparseArray[RawHamBloch/.rulef2];

Sort[Eigenvalues[ham2,3]]

{-4.2 21599673358501`,4.220699641899648`,4.221190634415838`}

Sort[Eigenvalues[ham2,10]]

{-4.221 681225493652`,-4.220934865650206`,-4.220621455741591`,-4.220023137110047`,-4.219207453227703`,4.21879283230247`,4.219605649533716`,4.220988985077544`,4.221062566397647`,4.221346747684675`}

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

{-4.2 2168122762485`,-4.220949631445061`,-4.220722846687062`,-4.220631517986182`,-4.220413040086286`,-4.220016414477176`,4.221017036779019`,4.221066732619922`,4.221281707170781`,4.221347789026139`}

liens[[1]]

Norm[cparmodcom[[1]]-cparmodcom[[11761]]]*locto[[7]]*Sqrt[2]

(*** Evaluation du temps de calcul des spectres complets en fonction de la taille du système ***)

rulef={Phi->1};
ham=SparseArray[RawHamFlux/.rulef];

eig3=Sort[Eigenvalues[ham]];