Hamiltonien sous champ magnétique, PBC sans phasons

Octo[1]=1;
Octo[2]=2;
Octo[n_]:=2Octo[n-1]+Octo[n-2];

Octo[4]

(*Importation des coordonnées des points et des liens*)

order=5; (*ordre de l'approximant*)
liens=Import["/Users/msoret/Desktop/Stage_M2/pavage_octo/data_pavage_octo/liensPBC"<>ToString[order]<>".dat"];
ppts=Import["/Users/msoret/Desktop/Stage_M2/pavage_octo/data_pavage_octo/PtsPBC"<>ToString[order]<>".dat"];
l=Length[liens];
ll=Length[ppts];
cut=0.51;

(*Construction du Hamiltonien*)

RawHamFlux=Reap[Do[{rat=Sqrt[2.]*Octo[order+1];
dx=ppts[[liens[[m,2]],1]]-ppts[[liens[[m,1]],1]];
dy=ppts[[liens[[m,2]],2]]-ppts[[liens[[m,1]],2]];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dx]<cut)&&(Abs[dy]<cut),(*Print["Cond  "," ",dx," ",dy," ",liens[[m,1]]," ",liens[[m,2]]];*)
x1=ppts[[liens[[m,1]],1]];
x2=ppts[[liens[[m,2]],1]];
y1=ppts[[liens[[m,1]],2]];
y2=ppts[[liens[[m,2]],2]];
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(dx>cut)&&(Abs[dy]<cut),(*1 proche du bord gauche, lien 1->2 croisant bord gauche*)
x1=ppts[[liens[[m,1]],1]];
x2=ppts[[liens[[m,2]],1]]-1;
y1=ppts[[liens[[m,1]],2]];
y2=ppts[[liens[[m,2]],2]];
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2+2*I*Pi*Phi*rat*rat*y2]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2-2*I*Pi*Phi*rat*rat*y2]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(dx<-cut)&&(Abs[dy]<cut),(*1 proche du bord droit, lien 1->2 croisant bord droit*)
x1=ppts[[liens[[m,1]],1]]-1;
x2=ppts[[liens[[m,2]],1]];
y1=ppts[[liens[[m,1]],2]];
y2=ppts[[liens[[m,2]],2]];
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2-2*Pi*Phi*I*rat*rat*y1]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2+2*Pi*Phi*I*rat*rat*y1]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dx]<cut)&&(dy>cut),(*1 proche du bord bas, lien 1->2 croisant bord bas*)
x1=ppts[[liens[[m,1]],1]];
x2=ppts[[liens[[m,2]],1]];
y1=ppts[[liens[[m,1]],2]];
y2=ppts[[liens[[m,2]],2]]-1;
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dx]<cut)&&(dy<-cut),(*1 proche du bord haut, lien 1->2 croisant bord haut*)
x1=ppts[[liens[[m,1]],1]];
x2=ppts[[liens[[m,2]],1]];
y1=ppts[[liens[[m,1]],2]]-1;
y2=ppts[[liens[[m,2]],2]];
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(dx>cut)&&(dy>cut),(*1 proche du coin bas gauche, lien 1->2 croisant bords gauche et bas*)
x1=ppts[[liens[[m,1]],1]];
x2=ppts[[liens[[m,2]],1]]-1;
y1=ppts[[liens[[m,1]],2]];
y2=ppts[[liens[[m,2]],2]]-1;
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2+2*I*Pi*Phi*rat*rat*(y2+1)]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2-2*I*Pi*Phi*rat*rat*(y2+1)]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(dx<-cut)&&(dy<-cut),(*1 proche du coin haut droit, lien 1->2 croisant bords haut et droit*)
x1=ppts[[liens[[m,1]],1]]-1;
x2=ppts[[liens[[m,2]],1]];
y1=ppts[[liens[[m,1]],2]]-1;
y2=ppts[[liens[[m,2]],2]];
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2-2*I*Pi*Phi*rat*rat*(y1+1)]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-2*I*Pi*Phi*rat*rat*(y2-y1)*(x2+x1)/2+2*I*Pi*Phi*rat*rat*(y1+1)]}];];

},{m,1,l}]][[2,1]];

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

Length[ppts]

Hamiltonien sous champ magnétique, OBC avec phasons

order=3; (*ordre de l'approximant*)
liens=Import["/Users/msoret/Desktop/Stage_M2/pavage_octo/data_pavage_octo/liensPBC"<>ToString[order]<>".dat"];
ppts=Import["/Users/msoret/Desktop/Stage_M2/pavage_octo/data_pavage_octo/PtsPBC"<>ToString[order]<>".dat"];
l=Length[liens];
ll=Length[ppts];
cut=0.51;

RawHamFlux=Reap[Do[{rat=N[Octo[order+1]*Sqrt[2.]];
dx=ppts[[liens[[m,1]],1]]-ppts[[liens[[m,2]],1]];
dy=ppts[[liens[[m,1]],2]]-ppts[[liens[[m,2]],2]];

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

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

Papillon sous champ magnétique, PBC

(* Sans phasons *)

spfluxpbc=Monitor[Reap[
For[n=1,n≤2*Octo[order+1]^2+1,n++,
rulef={Phi->(n-1)/(2*Octo[order+1]^2)};
ham=SparseArray[RawHamFlux/.rulef];
Sow[Flatten[{(n-1)/(2*Octo[order+1]^2),Sort[Eigenvalues[ham]]}]/.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];

pbc=Table[ListPlot[Transpose[{Transpose[spfluxpbc][[1]],Transpose[spfluxpbc][[i]]}],PlotStyle->{Blue,PointSize[(10^(-8))]},PlotRange->{{0,1},{-4.23,4.23}},Joined->False,Axes->False],{i,2,ll+1}];

papillon=Show[pbc,Frame->True,FrameStyle->Thickness[0.002]]

(* Avec phasons *)

spfluxphason=Monitor[Reap[
For[n=1,n≤2*Octo[k+1]^2+1,n++,
rulef={Phi->(n-1)/(2*Octo[k+1]^2+1)};
ham=SparseArray[RawHamFluxPhas/.rulef];
Sow[Flatten[{(n-1)/(2*Octo[k+1]^2+1),Sort[Eigenvalues[ham]]}]/.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];

pph=Table[ListPlot[Transpose[{Transpose[spfluxphason][[1]],Transpose[spfluxphason][[i]]}],PlotStyle->{Blue,PointSize[(10^(-8))]},PlotRange->{{0,1},{-4.23,4.23}},Joined->False,Axes->False],{i,2,ll+1}];

papillonphas=Show[pph,Frame->True,FrameStyle->Thickness[0.0002]]