

%---
f[x_,y_]=f0+fx x + fy y + (fxx x^2 + fyy y^2 + 2 fxy x y)/2
ax=1
ay=0
bx=1/2
by=Sqrt[3]/2
f0=f[0,0]
f1=f[ax,0]
f2=f[bx,by]
f3=f[-bx,by]
f4=f[-ax,0]
f5=f[-bx,-by]
f6=f[bx,-by]
Solve[{f6==%12, f5==%11, f4==%10, f3==%9, f2==%8, f1==%7},{f0,fx,fy,fxx,fyy,fxy},
{f6,f5,f4,f3,f2,f1}]
Solve[{f6==f[bx,-by], f5==f[-bx,-by], f4==f[-ax,0], f3==f[-bx,by], f2==f[bx,by]},{fx,fy,fxx,fyy,fxy},
{f6,f5,f4,f3,f2,f1}]
(2/3)(f1+f2+f3+f4+f5+f6-6 f0)
(1/2)(f1-f4)
(1/2)((f2-f3)+(f6-f5))
1/(2 Sqrt[3]) (f2+f3-f5-f6)

%---
% 1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6
dx2y2m = (2/3){
{-6,1,0,1,1,0,0,1,0,0,0,0,1,1,0,0},
{1,-6,1,0,1,1,0,0,0,0,0,0,0,1,1,0},
{0,1,-6,1,0,1,1,0,0,0,0,0,0,0,1,1},
{1,0,1,-6,0,0,1,1,0,0,0,0,1,0,0,1},
{1,1,0,0,-6,1,0,1,1,0,0,1,0,0,0,0},
{0,1,1,0,1,-6,1,0,1,1,0,0,0,0,0,0},
{0,0,1,1,0,1,-6,1,0,1,1,0,0,0,0,0},
{1,0,0,1,1,0,1,-6,0,0,1,1,0,0,0,0},
{0,0,0,0,1,1,0,0,-6,1,0,1,1,0,0,1},
{0,0,0,0,0,1,1,0,1,-6,1,0,1,1,0,0},
{0,0,0,0,0,0,1,1,0,1,-6,1,0,1,1,0},
{0,0,0,0,1,0,0,1,1,0,1,-6,0,0,1,1},
{1,0,0,1,0,0,0,0,1,1,0,0,-6,1,0,1},
{1,1,0,0,0,0,0,0,0,1,1,0,1,-6,1,0},
{0,1,1,0,0,0,0,0,0,0,1,1,0,1,-6,1},
{0,0,1,1,0,0,0,0,1,0,0,1,1,0,1,-6}}

% 1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6
dxm=(1/2){
{0,1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0},
{-1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0,-1,0,1,0,0,0,0,0,0,0,0,0,0,0,0},
{1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0,0,0,0,0,1,0,-1,0,0,0,0,0,0,0,0},
{0,0,0,0,-1,0,1,0,0,0,0,0,0,0,0,0},
{0,0,0,0,0,-1,0,1,0,0,0,0,0,0,0,0},
{0,0,0,0,1,0,-1,0,0,0,0,0,0,0,0,0},
{0,0,0,0,0,0,0,0,0,1,0,-1,0,0,0,0},
{0,0,0,0,0,0,0,0,-1,0,1,0,0,0,0,0},
{0,0,0,0,0,0,0,0,0,-1,0,1,0,0,0,0},
{0,0,0,0,0,0,0,0,1,0,-1,0,0,0,0,0},
{0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,-1},
{0,0,0,0,0,0,0,0,0,0,0,0,-1,0,1,0},
{0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,1},
{0,0,0,0,0,0,0,0,0,0,0,0,1,0,-1,0}}

% 1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6
dym=1/(2 Sqrt[3]) {
{0,0,0,0,1,0,0,1,0,0,0,0,-1,-1,0,0},
{0,0,0,0,1,1,0,0,0,0,0,0,0,-1,-1,0},
{0,0,0,0,0,1,1,0,0,0,0,0,0,0,-1,-1},
{0,0,0,0,0,0,1,1,0,0,0,0,-1,0,0,-1},
{-1,-1,0,0,0,0,0,0,1,0,0,1,0,0,0,0},
{0,-1,-1,0,0,0,0,0,1,1,0,0,0,0,0,0},
{0,0,-1,-1,0,0,0,0,0,1,1,0,0,0,0,0},
{-1,0,0,-1,0,0,0,0,0,0,1,1,0,0,0,0},
{0,0,0,0,-1,-1,0,0,0,0,0,0,1,0,0,1},
{0,0,0,0,0,-1,-1,0,0,0,0,0,1,1,0,0},
{0,0,0,0,0,0,-1,-1,0,0,0,0,0,1,1,0},
{0,0,0,0,-1,0,0,-1,0,0,0,0,0,0,1,1},
{1,0,0,1,0,0,0,0,-1,-1,0,0,0,0,0,0},
{1,1,0,0,0,0,0,0,0,-1,-1,0,0,0,0,0},
{0,1,1,0,0,0,0,0,0,0,-1,-1,0,0,0,0},
{0,0,1,1,0,0,0,0,-1,0,0,-1,0,0,0,0}}

ax=1
ay=0
bx=1/2
by=Sqrt[3]/2

Do[Do[x[j 4 + i+1]=(ax/4)i+(bx/4)j-3/4; y[j 4 + i+1]=(ay/4)i+(by/4)j-by/2,{j,0,3,1}],{i,0,3,1}]

v1[x_,y_]=128(x^2+y^2-.25)
v[x_,y_]=Min[{0,v1[x,y],v1[x-ax,y-ay],v1[x+ax,y+ay],v1[x+bx,y+by],v1[x-bx,y-by]}]
vm=DiagonalMatrix[Table[v[x[i],y[i]],{i,16}]]

h[kx_,ky_]=((kx^2+ky^2)IdentityMatrix[8]-dx2y2m/.25^2)/2 -I kx dxm/.25 -I ky dym/.25 + vm

Eigenvalues[N[h[0,0]]]
Out[268]= {-18.8495, -0.2207, -0.2207, 2.38685, 7.87722, 7.87722, 13.3333, 
 >    18.6667, 18.6667, 28.1889, 28.1889, 32.2207, 32.2207, 38.6005, 38.6005, 
 >    40.4626}

{vals,vecs}=Eigensystem[N[h[0,0]]]
sortedvals=Sort[vals]
sortedvecs={vecs[[Position[vals,sortedvals[[1]]][[1,1]]]],vecs[[Position[vals,sortedvals[[2]]][[1,1]]]],
vecs[[Position[vals,sortedvals[[3]]][[1,1]]]],vecs[[Position[vals,sortedvals[[4]]][[1,1]]]]}
roundedprob=N[Round[1000*Abs[sortedvecs]^2/Max[Abs[sortedvecs]^2]]/(1000 )]

e1band[kx_,ky_]:=Min[Re[Eigenvalues[N[h[kx,ky]]]]]
ContourPlot[e1band[kx,ky],{kx,-3,3},{ky,-3,3}]

e2band[kx_,ky_]:=Sort[Re[Eigenvalues[N[h[kx,ky]]]]][[2]]
ContourPlot[e2band[kx,ky],{kx,-3,3},{ky,-3,3}]

e3band[kx_,ky_]:=Sort[Re[Eigenvalues[N[h[kx,ky]]]]][[3]]
ContourPlot[e3band[kx,ky],{kx,-3,3},{ky,-3,3}]

e4band[kx_,ky_]:=Sort[Re[Eigenvalues[N[h[kx,ky]]]]][[4]]
ContourPlot[e4band[kx,ky],{kx,-3,3},{ky,-3,3}]

Plot[{e1band[0,kx],e2band[0,kx],e3band[0,kx],e4band[0,kx]},{kx,0,Pi}]

%--

junk[i_,j_]:={i+8 j +1,Mod[i+1,8]+8 j +1,Mod[i-1,8]+8 j +1,i+8 Mod[j+1,8] +1,
Mod[i-1,8]+8 Mod[j+1,8] +1,i+8 Mod[j-1,8] +1,Mod[i+1,8]+8 Mod[j-1,8] +1}

dx2y2=DiagonalMatrix[Table[-6,{64}]]
Do[Do[dx2y2[[i+8 j +1,Mod[i+1,8]+8 j +1]]=1;
dx2y2[[i+8 j +1,Mod[i-1,8]+8 j +1]]=1;
dx2y2[[i+8 j +1,i+8 Mod[j+1,8] +1]]=1;
dx2y2[[i+8 j +1,Mod[i-1,8]+8 Mod[j+1,8] +1]]=1;
dx2y2[[i+8 j +1,i+8 Mod[j-1,8] +1]]=1;
dx2y2[[i+8 j +1,Mod[i+1,8]+8 Mod[j-1,8] +1]]=1,{i,0,7,1}],{j,0,7,1}]

dx=DiagonalMatrix[Table[0,{64}]]
Do[Do[dx[[i+8 j +1,Mod[i+1,8]+8 j +1]]=1;
dx[[i+8 j +1,Mod[i-1,8]+8 j +1]]=-1,{i,0,7,1}],{j,0,7,1}]

dy=DiagonalMatrix[Table[0,{64}]]
Do[Do[dy[[i+8 j +1,i+8 Mod[j+1,8] +1]]=1;
dy[[i+8 j +1,Mod[i-1,8]+8 Mod[j+1,8] +1]]=1;
dy[[i+8 j +1,i+8 Mod[j-1,8] +1]]=-1;
dy[[i+8 j +1,Mod[i+1,8]+8 Mod[j-1,8] +1]]=-1,{i,0,7,1}],{j,0,7,1}]

ax=1
ay=0
bx=1/2
by=Sqrt[3]/2

Do[Do[x[j 8 + i+1]=(ax/8)i+(bx/8)j-3/4; y[j 8 + i+1]=(ay/8)i+(by/8)j-by/2,{j,0,7,1}],{i,0,7,1}]

v1[x_,y_]=256(x^2+y^2-.25)
v[x_,y_]=Min[{0,v1[x,y],v1[x-ax,y-ay],v1[x+ax,y+ay],v1[x+bx,y+by],v1[x-bx,y-by]}]
u=DiagonalMatrix[Table[v[x[i],y[i]],{i,64}]]

h[kx_,ky_]=( (kx^2+ky^2)IdentityMatrix[64]-(2/3)dx2y2/.125^2)/2 -I kx (1/2) dx/.125 -I ky 1/(2 Sqrt[3]) dy/.125 + u

{vals,vecs}=Eigensystem[N[h[0,1]]]
sortedvals=Sort[vals]
sortedvecs={vecs[[Position[vals,sortedvals[[1]]][[1,1]]]],vecs[[Position[vals,sortedvals[[2]]][[1,1]]]],
vecs[[Position[vals,sortedvals[[3]]][[1,1]]]],vecs[[Position[vals,sortedvals[[4]]][[1,1]]]]}
roundedprob=N[Round[1000*Abs[sortedvecs]^2/Max[Abs[sortedvecs]^2]]/(1000 )]



val=Table[0,{41}]
valx=Table[0,{41}]
valxy=Table[0,{21}]
Do[ val[[i+1]]=Sort[Re[Eigenvalues[N[h[0,(i/40) 2 Pi/Sqrt[3]]]]]],{i,0,40,1}]
Do[ valx[[i+1]]=Sort[Re[Eigenvalues[N[h[(i/40) Pi/(Sqrt[3]/2)^2,0]]]]],{i,0,40,1}]
Do[ valxy[[i+1]]=Sort[Re[Eigenvalues[N[h[(i/20) (1/2) Pi/(Sqrt[3]/2)^2,2 Pi/Sqrt[3]]]]]],{i,0,20,1}]

<<Graphics`MultipleListPlot`
MultipleListPlot[Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[1]]},{i,41}],Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[2]]},{i,41}],
Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[3]]},{i,41}],Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[4]]},{i,41}],
Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[5]]},{i,41}],Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[6]]},{i,41}],
Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[7]]},{i,41}],
PlotJoined->True]
Show[
Plot[{ekfree[0,ky,7],ekfree[0,ky,1],ekfree[0,ky,2],ekfree[0,ky,3],ekfree[0,ky,4],ekfree[0,ky,5],ekfree[0,ky,6]},
{ky,0,2 Pi/Sqrt[3]},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]}}],
Plot[ekfree[0,ky,8],{ky,2,2 Pi/Sqrt[3]},PlotStyle->{{RGBColor[0,0,1]}}],
ListPlot[Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[1]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[2]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[3]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[4]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[5]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[6]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[7]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}]
]

MultipleListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[1]]},{i,41}],Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[2]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[3]]},{i,41}],Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[4]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[5]]},{i,41}],Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[6]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[7]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[8]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[9]]},{i,41}],
PlotJoined->True]
Show[%,Plot[ekfree[kx,0,10],{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[0,0,1]}}],
Plot[ekfree[kx,0,9],{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[1,0,1]}}],
Plot[ekfree[kx,0,8],{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[0,0,1]}}]]

Show[
Plot[{ekfree[kx,0,7],ekfree[kx,0,1],ekfree[kx,0,2],ekfree[kx,0,3],ekfree[kx,0,4],ekfree[kx,0,5],ekfree[kx,0,6]},
{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]}}],
Plot[ekfree[kx,0,10],{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[0,0,1]}}],
Plot[ekfree[kx,0,9],{kx,3,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[1,0,1]}}],
Plot[ekfree[kx,0,8],{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[0,0,1]}}],
ListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[1]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[2]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[3]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[4]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[5]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[6]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[7]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[8]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[9]]},{i,41}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}]
]

MultipleListPlot[Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[1]]},{i,0,20,1}],Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[2]]},{i,0,20,1}],
Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[3]]},{i,0,20,1}],Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[4]]},{i,0,20,1}],
Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[5]]},{i,0,20,1}],Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[6]]},{i,0,20,1}],
Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[7]]},{i,0,20,1}],
PlotJoined->True]

Show[
Plot[{ekfree[kx,2 Pi/Sqrt[3],2],ekfree[kx,2 Pi/Sqrt[3],3],
ekfree[kx,2 Pi/Sqrt[3],4],ekfree[kx,2 Pi/Sqrt[3],5],ekfree[kx,2 Pi/Sqrt[3],6],ekfree[kx,2 Pi/Sqrt[3],7]},
{kx,0,(1/2)Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[1,0,0]}}],
Plot[{ekfree[kx,2 Pi/Sqrt[3],8],ekfree[kx,2 Pi/Sqrt[3],11]},{kx,0,(1/2)Pi/(Sqrt[3]/2)^2},
PlotStyle->{{RGBColor[0,0,1]},{RGBColor[0,0,1]}}],
ListPlot[Table[{(i/20) (1/2) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[1]]},{i,0,20,1}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{(i/20) (1/2) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[2]]},{i,0,20,1}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{(i/20) (1/2) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[3]]},{i,0,20,1}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{(i/20) (1/2) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[4]]},{i,0,20,1}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{(i/20) (1/2) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[5]]},{i,0,20,1}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{(i/20) (1/2) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[6]]},{i,0,20,1}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}],
ListPlot[Table[{(i/20) (1/2) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[7]]},{i,0,20,1}],PlotJoined->True,PlotStyle->{AbsolutePointSize[0]}]]




v1[x_,y_]=128(x^2+y^2-.25)
v[x_,y_]=Min[{0,v1[x,y],v1[x-ax,y-ay],v1[x+ax,y+ay],v1[x+bx,y+by],v1[x-bx,y-by]}]
u=DiagonalMatrix[Table[v[x[i],y[i]],{i,64}]]
h[kx_,ky_]=( (kx^2+ky^2)IdentityMatrix[64]-(2/3)dx2y2/.125^2)/2 -I kx (1/2) dx/.125 -I ky 1/(2 Sqrt[3]) dy/.125 + u

Table[{(i/40) 2 Pi/Sqrt[3],val[[i]][[1]]/((i/40) 2 Pi/Sqrt[3])^2},{i,1,9,1}]
Table[{(i/40) 2 Pi/Sqrt[3],valx[[i]][[1]]/((i/40) 2 Pi/Sqrt[3])^2},{i,1,9,1}]

Table[{(i/40) 2 Pi/Sqrt[3],If[i<9,val[[i]][[2]],val[[i]][[1]]]},{i,41}]
Fit[%,{1,x,x^2},x]
Table[{(i/40) 2 Pi/Sqrt[3],If[i<9,valx[[i]][[2]],valx[[i]][[1]]]},{i,41}]
Fit[%,{1,x,x^2},x]

Table[{(i/40) 2 Pi/Sqrt[3],If[i<10,val[[i]][[4]],valx[[i]][[3]]]},{i,41}]
Fit[%,{1,x,x^2},x]

kax=0; kay=2 Pi/(Sqrt[3]/2)
kbx=kay Cos[Pi/6]; kby=-kay Sin[Pi/6]

crot=Exp[-I Pi/3]
klattice=Table[{0,0},{11}]
Do[klattice[[i]]={Re[2 Pi/(Sqrt[3]/2) I crot^(i-1)],Im[2 Pi/(Sqrt[3]/2) I crot^(i-1)]},{i,1,6}]
klattice[[8]]=-(klattice[[1]]+klattice[[2]])
klattice[[9]]=-(klattice[[2]]+klattice[[2]])
klattice[[10]]=-(klattice[[3]]+klattice[[2]])
klattice[[11]]=-(klattice[[1]]+klattice[[6]])

ekfree[kx_,ky_,i_]:=(1/2)( (kx+klattice[[i]][[1]])^2 +(ky+klattice[[i]][[2]])^2 )

%exact solution to free discrete matrix:

hex[h_]=2Cos[kx h]+4 Cos[kx h/2]Cos[ky Sqrt[3] h/2]-6
Series[hex[x],{x,0,6}]
Simplify[%]

              2     2   2        2     2 2  4
        -3 (kx  + ky ) x    3 (kx  + ky )  x
Out[4]= ----------------- + ----------------- + 
                2                  32
 
            6        4   2        2   4       6   6
     (-11 kx  - 15 kx  ky  - 45 kx  ky  - 9 ky ) x        7
>    ---------------------------------------------- + O[x]
                          3840



hexact8[kx_,ky_,i_]:=( (kx^2+ky^2)-(2/3)(2 Cos[klattice[[i]][[1]] .125] + 
4 Cos[klattice[[i]][[1]] .125/2] Cos[klattice[[i]][[2]] Sqrt[3] .125/2]-6)/.125^2)/2 + 
kx (1/2) (2 Sin[klattice[[i]][[1]] .125])/.125 + 
ky 1/(2 Sqrt[3]) (2 Cos[klattice[[i]][[1]] .125/2] 2 Sin[klattice[[i]][[2]] Sqrt[3] .125/2])/.125 

Plot[{hexact8[x 4 Pi/3,0,10],hexact8[x 2 Pi, x 2 Pi/Sqrt[3],8]},{x,0,1}]
Show[
Plot[{ekfree[kx,0,7],ekfree[kx,0,1],ekfree[kx,0,2],ekfree[kx,0,3],ekfree[kx,0,4],ekfree[kx,0,5],ekfree[kx,0,6]},
{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]}}],
Plot[ekfree[kx,0,10],{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[0,0,1]}}],
Plot[ekfree[kx,0,9],{kx,3,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[1,0,1]}}],
Plot[ekfree[kx,0,8],{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[0,0,1]}}],
Plot[{hexact8[kx,0,7],hexact8[kx,0,1],hexact8[kx,0,2],hexact8[kx,0,3],hexact8[kx,0,4],hexact8[kx,0,5],hexact8[kx,0,6]},
{kx,0,Pi/(Sqrt[3]/2)^2}],
Plot[hexact8[kx,0,10],{kx,0,Pi/(Sqrt[3]/2)^2}],
Plot[hexact8[kx,0,9],{kx,3,Pi/(Sqrt[3]/2)^2}],
Plot[hexact8[kx,0,8],{kx,0,Pi/(Sqrt[3]/2)^2}]
]

hexact16[kx_,ky_,i_]:=( (kx^2+ky^2)-(2/3)(2 Cos[klattice[[i]][[1]] (1/16)] + 
4 Cos[klattice[[i]][[1]] (1/16)/2] Cos[klattice[[i]][[2]] Sqrt[3] (1/16)/2]-6)/(1/16)^2)/2 + 
kx (1/2) (2 Sin[klattice[[i]][[1]] (1/16)])/(1/16) + 
ky 1/(2 Sqrt[3]) (2 Cos[klattice[[i]][[1]] (1/16)/2] 2 Sin[klattice[[i]][[2]] Sqrt[3] (1/16)/2])/(1/16) 

Show[
Plot[{ekfree[kx,0,7],ekfree[kx,0,1],ekfree[kx,0,2],ekfree[kx,0,3],ekfree[kx,0,4],ekfree[kx,0,5],ekfree[kx,0,6]},
{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]}}],
Plot[ekfree[kx,0,10],{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[0,0,1]}}],
Plot[ekfree[kx,0,9],{kx,3,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[1,0,1]}}],
Plot[ekfree[kx,0,8],{kx,0,Pi/(Sqrt[3]/2)^2},PlotStyle->{{RGBColor[0,0,1]}}],
Plot[{hexact16[kx,0,7],hexact16[kx,0,1],hexact16[kx,0,2],hexact16[kx,0,3],hexact16[kx,0,4],hexact16[kx,0,5],hexact16[kx,0,6]},
{kx,0,Pi/(Sqrt[3]/2)^2}],
Plot[hexact16[kx,0,10],{kx,0,Pi/(Sqrt[3]/2)^2}],
Plot[hexact16[kx,0,9],{kx,3,Pi/(Sqrt[3]/2)^2}],
Plot[hexact16[kx,0,8],{kx,0,Pi/(Sqrt[3]/2)^2}]
]

Do[Print[N[(Sum[Part[(klattice[[i]]+{0,2 Pi/Sqrt[3]})^(2),j],{j,2}])/2]],{i,10}]

Do[k0=(i/40) 2 Pi/Sqrt[3]; Print[N[k0,3]," ", N[ekfree[0,k0,7],3]," ", N[ekfree[0,k0,1],3]," ",
N[ekfree[0,k0,2],3]," ", N[ekfree[0,k0,3],3]," ", N[ekfree[0,k0,4],3]," ", N[ekfree[0,k0,8],3]],{i,0,40,1}]

Do[k0=(i/40) 4 Pi/3; Print[N[k0,3]," ", N[ekfree[k0,0,7],3]," ", N[ekfree[k0,0,2],3]," ",
N[ekfree[k0,0,1],3]," ", N[ekfree[k0,0,6],3]," ", N[ekfree[k0,0,10],3]],{i,0,40,1}]

Do[k0=(i/40) 2 Pi/3; Print[N[k0,3]," ", N[ekfree[k0,2 Pi/Sqrt[3],7],3]," ", N[ekfree[k0,2 Pi/Sqrt[3],6],3]," ",
N[ekfree[k0,2 Pi/Sqrt[3],3],3]," ", N[ekfree[k0,2 Pi/Sqrt[3],5],3]," ", N[ekfree[k0,2 Pi/Sqrt[3],2],3],
" ", N[ekfree[k0,2 Pi/Sqrt[3],1],3]],{i,0,40,1}]


Plot3D[ekfree[x,y,7],{x,2,4 Pi/3},{y,-1,1},ViewPoint->{-1,0,0}]
Plot3D[ekfree[x,y,5],{x,2,4 Pi/3},{y,-1,1},ViewPoint->{-1,0,0}]
Plot3D[ekfree[x,y,6],{x,2,4 Pi/3},{y,-1,1},ViewPoint->{-1,0,0}]
Show[%,%%,%%%,LightSources->{{{0,2,2},RGBColor[1,0,0]},{{+2,2,-2},RGBColor[0,1,0]},{{-2,2,-2},RGBColor[1,0,0]}}]
Show[%56,LightSources->{{{0,2,2},RGBColor[1,0,0]},{{+2,-2,2},RGBColor[0,1,0]},{{-2,-2,2},RGBColor[1,0,0]}}]

---
exact[kx_,ky_,gx_,g_y,hh_]:=( (kx^2+ky^2)-(2/3)(2 Cos[gx hh] + 
4 Cos[gx hh/2] Cos[gy Sqrt[3] hh/2]-6)/hh^2)/2 + 
kx (1/2) (2 Sin[gx hh])/hh + 
ky 1/(2 Sqrt[3]) (2 Cos[gx hh/2] 2 Sin[gy Sqrt[3] hh/2])/hh 


---

v1[x_,y_]=256(x^2+y^2-.25)
v[x_,y_]=Min[{0,v1[x,y],v1[x-ax,y-ay],v1[x+ax,y+ay],v1[x+bx,y+by],v1[x-bx,y-by]}]
u=DiagonalMatrix[Table[v[x[i],y[i]],{i,64}]]

h[kx_,ky_]=( (kx^2+ky^2)IdentityMatrix[64]-(2/3)dx2y2/.125^2)/2 -I kx (1/2) dx/.125 -I ky 1/(2 Sqrt[3]) dy/.125 + u

val=Table[0,{41}]
valx=Table[0,{41}]
valxy=Table[0,{21}]
Do[ val[[i+1]]=Sort[Re[Eigenvalues[N[h[0,(i/40) 2 Pi/Sqrt[3]]]]]],{i,0,40,1}]
Do[ valx[[i+1]]=Sort[Re[Eigenvalues[N[h[(i/40) Pi/(Sqrt[3]/2)^2,0]]]]],{i,0,40,1}]
Do[ valxy[[i+1]]=Sort[Re[Eigenvalues[N[h[(i/20) (1/2) Pi/(Sqrt[3]/2)^2,2 Pi/Sqrt[3]]]]]],{i,0,20,1}]

Do[k0=((i-1)/40) 2 Pi/Sqrt[3]; Print[N[k0,3]," ", N[val[[i]][[1]],3]," ", N[val[[i]][[2]],3]," ",
N[val[[i]][[3]],3]," ", N[val[[i]][[4]],3]," ", N[val[[i]][[5]],3]," ", N[val[[i]][[6]],3]," ",
N[val[[i]][[7]],3]],{i,41}]

Do[k0=((i-1)/40) 4 Pi/3; Print[N[k0,3]," ", N[valx[[i]][[1]],3]," ", N[valx[[i]][[2]],3]," ",
N[valx[[i]][[3]],3]," ", N[valx[[i]][[4]],3]," ", N[valx[[i]][[5]],3]," ", N[valx[[i]][[6]],3]," ",
 N[valx[[i]][[7]],3]],{i,41}]

Do[k0=((i-1)/20) 2 Pi/3; Print[N[k0,3]," ", N[valxy[[i]][[1]],3]," ", N[valxy[[i]][[2]],3]," ",
N[valxy[[i]][[3]],3]," ", N[valxy[[i]][[4]],3]," ", N[valxy[[i]][[5]],3]," ", N[valxy[[i]][[6]],3] ],{i,21}]

----

dx2y2=DiagonalMatrix[Table[-6,{64 4}]]
Do[Do[dx2y2[[i+16 j +1,Mod[i+1,16]+16 j +1]]=1;
dx2y2[[i+16 j +1,Mod[i-1,16]+16 j +1]]=1;
dx2y2[[i+16 j +1,i+16 Mod[j+1,16] +1]]=1;
dx2y2[[i+16 j +1,Mod[i-1,16]+16 Mod[j+1,16] +1]]=1;
dx2y2[[i+16 j +1,i+16 Mod[j-1,16] +1]]=1;
dx2y2[[i+16 j +1,Mod[i+1,16]+16 Mod[j-1,16] +1]]=1,{i,0,15,1}],{j,0,15,1}]

dx=DiagonalMatrix[Table[0,{64 4}]]
Do[Do[dx[[i+16 j +1,Mod[i+1,16]+16 j +1]]=1;
dx[[i+16 j +1,Mod[i-1,16]+16 j +1]]=-1,{i,0,15,1}],{j,0,15,1}]

dy=DiagonalMatrix[Table[0,{64 4}]]
Do[Do[dy[[i+16 j +1,i+16 Mod[j+1,16] +1]]=1;
dy[[i+16 j +1,Mod[i-1,16]+16 Mod[j+1,16] +1]]=1;
dy[[i+16 j +1,i+16 Mod[j-1,16] +1]]=-1;
dy[[i+16 j +1,Mod[i+1,16]+16 Mod[j-1,16] +1]]=-1,{i,0,15,1}],{j,0,15,1}]

ax=1
ay=0
bx=1/2
by=Sqrt[3]/2

Do[Do[x[j 16 + i+1]=(ax/16)i+(bx/16)j-3/4; y[j 16 + i+1]=(ay/16)i+(by/16)j-by/2,{j,0,15,1}],{i,0,15,1}]

v1[x_,y_]=256(x^2+y^2-.25)
v[x_,y_]=Min[{0,v1[x,y],v1[x-ax,y-ay],v1[x+ax,y+ay],v1[x+bx,y+by],v1[x-bx,y-by]}]
u=DiagonalMatrix[Table[v[x[i],y[i]],{i,64 4}]]

h[kx_,ky_]=( (kx^2+ky^2)IdentityMatrix[64 4]-(2/3)dx2y2/(1/16)^2)/2 -I kx (1/2) dx/(1/16) -I ky 1/(2 Sqrt[3]) dy/(1/16) + u

val=Table[0,{41}]
valx=Table[0,{41}]
valxy=Table[0,{21}]
Do[ val[[i+1]]=Sort[Re[Eigenvalues[N[h[0,(i/40) 2 Pi/Sqrt[3]]]]]]; Print[i],{i,0,40,1}]
Do[ valx[[i+1]]=Sort[Re[Eigenvalues[N[h[(i/40) Pi/(Sqrt[3]/2)^2,0]]]]]; Print[i],{i,0,40,1}]
Do[ valxy[[i+1]]=Sort[Re[Eigenvalues[N[h[(i/20) (1/2) Pi/(Sqrt[3]/2)^2,2 Pi/Sqrt[3]]]]]]; Print[i],{i,0,20,1}]

---
v1[x_,y_]=64(x^2+y^2-.25)
v[x_,y_]=Min[{0,v1[x,y],v1[x-ax,y-ay],v1[x+ax,y+ay],v1[x+bx,y+by],v1[x-bx,y-by]}]
u=DiagonalMatrix[Table[v[x[i],y[i]],{i,64}]]
h[kx_,ky_]=( (kx^2+ky^2)IdentityMatrix[64]-(2/3)dx2y2/.125^2)/2 -I kx (1/2) dx/.125 -I ky 1/(2 Sqrt[3]) dy/.125 + u
{vals,vecs}=Eigensystem[N[h[0,0]]]
sortedvals=Sort[vals]
sortedvecs={vecs[[Position[vals,sortedvals[[1]]][[1,1]]]],vecs[[Position[vals,sortedvals[[2]]][[1,1]]]],
vecs[[Position[vals,sortedvals[[3]]][[1,1]]]],vecs[[Position[vals,sortedvals[[4]]][[1,1]]]]}
rprob=N[Round[1000*Abs[sortedvecs]^2/Max[Abs[sortedvecs]^2]]/(1000 )]


Do[Print[rprob[[1]][[8 i+1]]," ",  rprob[[1]][[8 i+2]]," ",
rprob[[1]][[8 i+3]]," ", rprob[[1]][[8 i+4]]," ", rprob[[1]][[8 i+5]]," ", rprob[[1]][[8 i+6]]," ",
rprob[[1]][[8 i+7]]," ",rprob[[1]][[8 i+8]]],{i,0,7,1}]

Do[Print[rprob[[2]][[8 i+1]]," ",  rprob[[2]][[8 i+2]]," ",
rprob[[2]][[8 i+3]]," ", rprob[[2]][[8 i+4]]," ", rprob[[2]][[8 i+5]]," ", rprob[[2]][[8 i+6]]," ",
rprob[[2]][[8 i+7]]," ",rprob[[2]][[8 i+8]]],{i,0,7,1}]

Do[Print[rprob[[3]][[8 i+1]]," ",  rprob[[3]][[8 i+2]]," ",
rprob[[3]][[8 i+3]]," ", rprob[[3]][[8 i+4]]," ", rprob[[3]][[8 i+5]]," ", rprob[[3]][[8 i+6]]," ",
rprob[[3]][[8 i+7]]," ",rprob[[3]][[8 i+8]]],{i,0,7,1}]

Do[Print[rprob[[4]][[8 i+1]]," ", rprob[[4]][[8 i+2]]," ",
rprob[[4]][[8 i+3]]," ", rprob[[4]][[8 i+4]]," ", rprob[[4]][[8 i+5]]," ", rprob[[4]][[8 i+6]]," ",
rprob[[4]][[8 i+7]]," ",rprob[[4]][[8 i+8]]],{i,0,7,1}]

---
crot=Exp[-I Pi/3]
klattice=Table[{0,0},{11}]
Do[klattice[[i]]={Re[2 Pi/(Sqrt[3]/2) I crot^(i-1)],Im[2 Pi/(Sqrt[3]/2) I crot^(i-1)]},{i,1,6}]
klattice[[8]]=-(klattice[[1]]+klattice[[2]])
klattice[[9]]=-(klattice[[2]]+klattice[[2]])
klattice[[10]]=-(klattice[[3]]+klattice[[2]])
klattice[[11]]=-(klattice[[1]]+klattice[[6]])

u62p[x_,y_]=Exp[- I klattice[[6]].{x,y}]+Exp[- I klattice[[2]].{x,y}]
tvec=Table[u62p[x[i],y[i]],{i,64}]
trprob=N[Round[1000*Abs[tvec]^2/Max[Abs[tvec]^2]]/(1000 )]
Do[Print[trprob[[8 i+1]]," ", trprob[[8 i+2]]," ",
trprob[[8 i+3]]," ", trprob[[8 i+4]]," ", trprob[[8 i+5]]," ", trprob[[8 i+6]]," ",
trprob[[8 i+7]]," ",trprob[[8 i+8]]],{i,0,7,1}]

u62m[x_,y_]=Exp[- I klattice[[6]].{x,y}]-Exp[- I klattice[[2]].{x,y}]
tvec=Table[u62m[x[i],y[i]],{i,64}]

u71p[x_,y_]=1+Exp[- I klattice[[1]].{x,y}]
tvec=Table[u71p[x[i],y[i]],{i,64}]
u71m[x_,y_]=1-Exp[- I klattice[[1]].{x,y}]
tvec=Table[u71m[x[i],y[i]],{i,64}]

u32p[x_,y_]=Exp[- I klattice[[3]].{x,y}]+Exp[- I klattice[[2]].{x,y}]
tvec=Table[u32p[x[i],y[i]],{i,64}]
u32m[x_,y_]=Exp[- I klattice[[3]].{x,y}]-Exp[- I klattice[[2]].{x,y}]
tvec=Table[u32m[x[i],y[i]],{i,64}]
u41m[x_,y_]=Exp[- I klattice[[4]].{x,y}]-Exp[- I klattice[[1]].{x,y}]
tvec=Table[u41m[x[i],y[i]],{i,64}]

u732p[x_,y_]=1+Exp[- I klattice[[3]].{x,y}]+Exp[- I klattice[[2]].{x,y}]
tvec=Table[u732p[x[i],y[i]],{i,64}]
u732m[x_,y_]=-2+Exp[- I klattice[[3]].{x,y}]+Exp[- I klattice[[2]].{x,y}]
tvec=Table[u732m[x[i],y[i]],{i,64}]

uall[x_,y_]=Exp[- I klattice[[1]].{x,y}]+Exp[- I klattice[[2]].{x,y}]+
Exp[- I klattice[[3]].{x,y}]+Exp[- I klattice[[4]].{x,y}]+
Exp[- I klattice[[5]].{x,y}]+Exp[- I klattice[[6]].{x,y}]
tvec=Table[uall[x[i],y[i]],{i,64}]

q=.29816

ug[x_,y_]=(2+q)Exp[- I klattice[[1]].{x,y}]-Exp[- I klattice[[2]].{x,y}]-
(3+q)Exp[- I klattice[[3]].{x,y}]-(2+q)Exp[- I klattice[[4]].{x,y}]+
Exp[- I klattice[[5]].{x,y}]+(3+q)Exp[- I klattice[[6]].{x,y}]
tvec=Table[ug[x[i],y[i]],{i,64}]

q=.36039
ug2[x_,y_]=-(3+q)Exp[- I klattice[[1]].{x,y}]-(4+q)Exp[- I klattice[[2]].{x,y}]-
Exp[- I klattice[[3]].{x,y}]+(3+q)Exp[- I klattice[[4]].{x,y}]+
(4+q)Exp[- I klattice[[5]].{x,y}]+Exp[- I klattice[[6]].{x,y}]
tvec=Table[ug2[x[i],y[i]],{i,64}]

Table[Exp[I klattice[[6]].{x[i],y[i]}],{i,64}].Table[Exp[-I klattice[[1]].{x[i],y[i]}],{i,64}]

----
random point:{0.698132, 2.4184}
-----
{1, a, b, c, b, a}, {a, 1, a, b, c, b}, {b, a, 1, a, b, c}, 
 {c, b, a, 1, a, b}, {b, c, b, a, 1, a}, {a, b, c, b, a, 1}}

In[19]:= vals

Out[19]= {1 + a - b - c, 1 + a - b - c, 1 - 2 a + 2 b - c, 1 - a - b + c, 
 
>    1 - a - b + c, 1 + 2 a + 2 b + c}

In[20]:= vecs[[1]]

Out[20]= {1, 0, -1, -1, 0, 1}

In[21]:= vecs[[2]]

Out[21]= {-1, -1, 0, 1, 1, 0}

In[22]:= vecs[[4]]

Out[22]= {-1, 0, 1, -1, 0, 1}

In[23]:= vecs[[5]]

Out[23]= {-1, 1, 0, -1, 1, 0}

In[24]:= m.{1,1,1,1,1,1}

Solve[Det[{{x,y,z,t,u,v},{1, 0, -1, -1, 0, 1},{-1, -2, -1, 1, 2, 1},
{-1, 0, 1, -1, 0, 1},{-3, 2, 1, -3, 2, 1},{1,1,1,1,1,1}}]==0,{x,y,z,t,u,v}]
Solve[Det[{{-5,-1,1,-1,1,-1},{1, 0, -1, -1, 0, 1},{-1, -1, 0, 1, 1, 0},
{-1, 0, 1, -1, 0, 1},{-1, 1, 0, -1, 1, 0},{1,1,1,1,1,1}}]==0,{x,y,z,t,u,v}]

{-1, -2, -1, 1, 2, 1}
{-3, 2, 1, -3, 2, 1}



-----
Integrate[Exp[I k r Cos[t]],{t,-Pi,Pi}]
                                                2  2
Out[43]= If[Im[k r] == 0, 2 Pi BesselJ[0, Sqrt[k  r ]], 
 
                I k r Cos[t]
>    Integrate[E            , {t, -Pi, Pi}]]

Integrate[BesselJ[0, k r] (r^2-1/4)r,{r,0,1/2}]

                                     2                        2
                 2             Sqrt[k ]                 Sqrt[k ]
Out[45]= -(Sqrt[k ] BesselJ[1, --------] - 4 BesselJ[2, --------] + 
                                  2                        2
 
                                  2
              2             Sqrt[k ]         2
>       Sqrt[k ] BesselJ[3, --------]) / (8 k )
                               2
  - 2 Pi (k*BesselJ[1, k/2] - 4*BesselJ[2, k/2] + 
      k*BesselJ[3, k/2])/(8*k^2)


In[54]:=  Series[  - 2 Pi (k*BesselJ[1, k/2] - 4*BesselJ[2, k/2] + 
      k*BesselJ[3, k/2])/(8*k^2) ,{k,0,4}]

                   2       4
         -Pi   Pi k    Pi k         5
Out[54]= --- + ----- - ------ + O[k]
         32    1536    196608

Plot[- 2 Pi (k*BesselJ[1, k/2] - 4*BesselJ[2, k/2] + 
      k*BesselJ[3, k/2])/(8*k^2) ,{k,0,30},PlotRange->{-.1,.01}]

Integrate[2 Pi (r^2-1/4) r,{r,0,1/2}]

----
crot=Exp[-I Pi/3]
klattice=Table[{0,0},{37}]
Do[klattice[[i+1]]={Re[2 Pi/(Sqrt[3]/2) I crot^(i-1)],Im[2 Pi/(Sqrt[3]/2) I crot^(i-1)]},{i,1,6}]
temp=(klattice[[2]]+klattice[[3]])
Do[klattice[[i+7]]={Re[(temp[[1]]+I temp[[2]]) crot^(i-1)],Im[(temp[[1]]+I temp[[2]]) crot^(i-1)]},{i,1,6}]
temp=2 klattice[[2]]
Do[klattice[[i+13]]={Re[(temp[[1]]+I temp[[2]]) crot^(i-1)],Im[(temp[[1]]+I temp[[2]]) crot^(i-1)]},{i,1,6}]
temp=(2 klattice[[2]]+klattice[[3]])
Do[klattice[[(2 i-1)+19]]={Re[(temp[[1]]+I temp[[2]]) crot^(i-1)],Im[(temp[[1]]+I temp[[2]]) crot^(i-1)]},{i,1,6}]
temp=(klattice[[2]]+2 klattice[[3]])
Do[klattice[[(2 i)+19]]={Re[(temp[[1]]+I temp[[2]]) crot^(i-1)],Im[(temp[[1]]+I temp[[2]]) crot^(i-1)]},{i,1,6}]
temp=3 klattice[[2]]
Do[klattice[[i+31]]={Re[(temp[[1]]+I temp[[2]]) crot^(i-1)],Im[(temp[[1]]+I temp[[2]]) crot^(i-1)]},{i,1,6}]

ListPlot[klattice,AspectRatio->Automatic]

v0=1

vk[k_]=v0(- 2 Pi (k*BesselJ[1, k/2] - 4*BesselJ[2, k/2] +  k*BesselJ[3, k/2])/(8*k^2))/(Sqrt[3]/2)

vkm=DiagonalMatrix[Table[-v0 Pi/32 /(Sqrt[3]/2),{i,37}]]

Do[Do[tempv=vk[Sqrt[Sum[Part[(klattice[[i]]-klattice[[j]])^(2),l],{l,2}]]];
vkm[[i,j]]=tempv; vkm[[j,i]]=tempv, {j,i+1,37}],{i,1,37}]


Table[Sum[Part[(klattice[[i]]-klattice[[j]])^(2),l],{l,2}],{j,1,37},{i,1,37}]
Simplify[%]
Flatten[%]
Union[%]
klist=Sqrt[%]
ListPlot[Table[{klist[[l]],vk[klist[[l]]]},{l,15}]]
Plot[vk[k],{k,0,45}]
Show[%,%%,PlotRange->{-.031,.007}]

h[kx_,ky_]=N[(1/2)DiagonalMatrix[Table[(kx+klattice[[i]][[1]])^2+(ky+klattice[[i]][[2]])^2,{i,37}]]+vkm]

val=Table[0,{41}]
valx=Table[0,{41}]
valxy=Table[0,{21}]
Do[ val[[i+1]]=Sort[Re[Eigenvalues[h[0,(i/40) 2 Pi/Sqrt[3]]]]]; Print[i],{i,0,40,1}]
Do[ valx[[i+1]]=Sort[Re[Eigenvalues[h[(i/40) Pi/(Sqrt[3]/2)^2,0]]]],{i,0,40,1}]
Do[ valxy[[i+1]]=Sort[Re[Eigenvalues[h[(i/20) (1/2) Pi/(Sqrt[3]/2)^2,2 Pi/Sqrt[3]]]]],{i,0,20,1}]

<<Graphics`MultipleListPlot`
MultipleListPlot[Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[1]]},{i,41}],Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[2]]},{i,41}],
Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[3]]},{i,41}],Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[4]]},{i,41}],
Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[5]]},{i,41}],Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[6]]},{i,41}],
Table[{((i-1)/40) 2 Pi/Sqrt[3],val[[i]][[7]]},{i,41}],
PlotJoined->True]
MultipleListPlot[Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[1]]},{i,41}],Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[2]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[3]]},{i,41}],Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[4]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[5]]},{i,41}],Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[6]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[7]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[8]]},{i,41}],
Table[{((i-1)/40) Pi/(Sqrt[3]/2)^2,valx[[i]][[9]]},{i,41}],
PlotJoined->True]
MultipleListPlot[Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[1]]},{i,0,20,1}],Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[2]]},{i,0,20,1}],
Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[3]]},{i,0,20,1}],Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[4]]},{i,0,20,1}],
Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[5]]},{i,0,20,1}],Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[6]]},{i,0,20,1}],
Table[{(i/40) Pi/(Sqrt[3]/2)^2,valxy[[i+1]][[7]]},{i,0,20,1}],
PlotJoined->True]

Do[k0=((i-1)/40) 2 Pi/Sqrt[3]; Print[N[k0,3]," ", N[val[[i]][[1]],3]," ", N[val[[i]][[2]],3]," ",
N[val[[i]][[3]],3]," ", N[val[[i]][[4]],3]," ", N[val[[i]][[5]],3]," ", N[val[[i]][[6]],3]," ",
N[val[[i]][[7]],3]],{i,41}]

Do[k0=((i-1)/40) 4 Pi/3; Print[N[k0,3]," ", N[valx[[i]][[1]],3]," ", N[valx[[i]][[2]],3]," ",
N[valx[[i]][[3]],3]," ", N[valx[[i]][[4]],3]," ", N[valx[[i]][[5]],3]," ", N[valx[[i]][[6]],3]," ",
 N[valx[[i]][[7]],3]],{i,41}]

Do[k0=((i-1)/20) 2 Pi/3; Print[N[k0,3]," ", N[valxy[[i]][[1]],3]," ", N[valxy[[i]][[2]],3]," ",
N[valxy[[i]][[3]],3]," ", N[valxy[[i]][[4]],3]," ", N[valxy[[i]][[5]],3]," ", N[valxy[[i]][[6]],3] ],{i,21}]

