# How to improve the speed of this numerical integration which includes calculations of large matrix?

The integrand is Re[k^2*Sin[2*theta]*T], where k is constant, theta is a variable (x, y and fi are the other variables), and T includes calculation of large matrix of M1 and M2. For the codes of original version, I found it needs more than 240s to work out the numerical integration of G (x0 and y0 are given). How to improve the speed of the numerical integration? Many thanks!

PS: Thanks very much for the answer provided by Henrik Schumacher. I have revised the code as the answer's suggestion, and the revised version is at last of the code. Can the speed be further improved?

original version:
Clear["Global*"]
x0 = 0;
y0 = -50;
afa = Pi/180*0.2;
K0 = 29.5;
kc = 5.12317;
k = kc*Sqrt[5];
v1 = kc^2*4;
v2 = kc^2*2;
v3 = kc^2*1;
r0 = {0.142, 0, 0};
a1 = {0.213, 0.123, 0};
a2 = {0.213, -0.123, 0};
b1 = {14.75, 25.55, 0};
b2 = {14.75, -25.55, 0};
Kh[0] = {0, 0, 0};
Ks[0] = {0, 0, 0};
Kh[1] = {K0, 0, 0};
Ks[1] = {K0*Cos[afa], K0*Sin[afa], 0};
Kh[2] = {K0/2, K0*Sqrt[3]/2, 0};
Ks[2] = {K0*Cos[afa + Pi/3], K0*Sin[afa + Pi/3], 0};
Kh[3] = {-K0/2, K0*Sqrt[3]/2, 0};
Ks[3] = {K0*Cos[afa + 2*Pi/3], K0*Sin[afa + 2*Pi/3], 0};
Kh[4] = {-K0, 0, 0};
Ks[4] = {K0*Cos[afa + Pi], K0*Sin[afa + Pi], 0};
Kh[5] = {-K0/2, -K0*Sqrt[3]/2, 0};
Ks[5] = {K0*Cos[afa + 4*Pi/3], K0*Sin[afa + 4*Pi/3], 0};
Kh[6] = {K0/2, -K0*Sqrt[3]/2, 0};
Ks[6] = {K0*Cos[afa + 5*Pi/3], K0*Sin[afa + 5*Pi/3], 0};
integrand[x_?NumericQ, y_, theta_, fi_] :=
Module[{ks, r, p, r1, x1, y1, r2, u, tranB, rB, tranA, rA, M1, M2,
S1, S2, S, T, S11, S12, S13, S14, S15, S16, S17, S21, S22, S23,
S24, S25, S26, S27, S31, S32, S33, S34, S35, S36, S37, S41, S42,
S43, S44, S45, S46, S47, S51, S52, S53, S54, S55, S56, S57, S61,
S62, S63, S64, S65, S66, S67, S71, S72, S73, S74, S75, S76, S77},
ks = {k*Sin[theta]*Cos[fi], k*Sin[theta]*Sin[fi], k*Cos[theta]};
r = {x, y, 0};
p[i_] := k^2 - Total[(ks + Kh[i])^2];
r1 = {-y*afa, x*afa, 0};
x1 = Mod[r1 . b1/(2*Pi), 1];
y1 = Mod[r1 . b2/(2*Pi), 1];
r2 = x1*a1 + y1*a2;
u = First[SortBy[{r0 - r2, 2*r0 - r2}, Norm]];
tranB =
Min[Norm[r2], Norm[a1 - r2], Norm[a2 - r2], Norm[2*r0 - r2],
Norm[3*r0 - r2], Norm[a1 - r0 - r2], Norm[a2 - r0 - r2]];
rB = Sqrt[tranB^2 + 0.112];
tranA =
Min[Norm[r2], Norm[a1 - r2], Norm[a2 - r2], Norm[r0 - r2],
Norm[3*r0 - r2], Norm[a1 + r0 - r2], Norm[a2 + r0 - r2]];
rA = Sqrt[tranA^2 + 0.112];
M1 = {{0, v1*(1 + Exp[I*((Kh[0] - Kh[1]) . r0)]),
v1*(1 + Exp[I*((Kh[0] - Kh[2]) . r0)]),
v1*(1 + Exp[I*((Kh[0] - Kh[3]) . r0)]),
v1*(1 + Exp[I*((Kh[0] - Kh[4]) . r0)]),
v1*(1 + Exp[I*((Kh[0] - Kh[5]) . r0)]),
v1*(1 + Exp[I*((Kh[0] - Kh[6]) . r0)])}, {v1*(1 +
Exp[I*((Kh[1] - Kh[0]) . r0)]), p[1],
v1*(1 + Exp[I*((Kh[1] - Kh[2]) . r0)]),
v2*(1 + Exp[I*((Kh[1] - Kh[3]) . r0)]),
v3*(1 + Exp[I*((Kh[1] - Kh[4]) . r0)]),
v2*(1 + Exp[I*((Kh[1] - Kh[5]) . r0)]),
v1*(1 + Exp[I*((Kh[1] - Kh[6]) . r0)])}, {v1*(1 +
Exp[I*((Kh[2] - Kh[0]) . r0)]),
v1*(1 + Exp[I*((Kh[2] - Kh[1]) . r0)]), p[2],
v1*(1 + Exp[I*((Kh[2] - Kh[3]) . r0)]),
v2*(1 + Exp[I*((Kh[2] - Kh[4]) . r0)]),
v3*(1 + Exp[I*((Kh[2] - Kh[5]) . r0)]),
v2*(1 + Exp[I*((Kh[2] - Kh[6]) . r0)])}, {v1*(1 +
Exp[I*((Kh[3] - Kh[0]) . r0)]),
v2*(1 + Exp[I*((Kh[3] - Kh[1]) . r0)]),
v1*(1 + Exp[I*((Kh[3] - Kh[2]) . r0)]), p[3],
v1*(1 + Exp[I*((Kh[3] - Kh[4]) . r0)]),
v2*(1 + Exp[I*((Kh[3] - Kh[5]) . r0)]),
v3*(1 + Exp[I*((Kh[3] - Kh[6]) . r0)])}, {v1*(1 +
Exp[I*((Kh[4] - Kh[0]) . r0)]),
v3*(1 + Exp[I*((Kh[4] - Kh[1]) . r0)]),
v2*(1 + Exp[I*((Kh[4] - Kh[2]) . r0)]),
v1*(1 + Exp[I*((Kh[4] - Kh[3]) . r0)]), p[4],
v1*(1 + Exp[I*((Kh[4] - Kh[5]) . r0)]),
v2*(1 + Exp[I*((Kh[4] - Kh[6]) . r0)])}, {v1*(1 +
Exp[I*((Kh[5] - Kh[0]) . r0)]),
v2*(1 + Exp[I*((Kh[5] - Kh[1]) . r0)]),
v3*(1 + Exp[I*((Kh[5] - Kh[2]) . r0)]),
v2*(1 + Exp[I*((Kh[5] - Kh[3]) . r0)]),
v1*(1 + Exp[I*((Kh[5] - Kh[4]) . r0)]), p[5],
v1*(1 + Exp[I*((Kh[5] - Kh[6]) . r0)])}, {v1*(1 +
Exp[I*((Kh[6] - Kh[0]) . r0)]),
v1*(1 + Exp[I*((Kh[6] - Kh[1]) . r0)]),
v2*(1 + Exp[I*((Kh[6] - Kh[2]) . r0)]),
v3*(1 + Exp[I*((Kh[6] - Kh[3]) . r0)]),
v2*(1 + Exp[I*((Kh[6] - Kh[4]) . r0)]),
v1*(1 + Exp[I*((Kh[6] - Kh[5]) . r0)]), p[6]}};
M2 = {{0,
v1*(Exp[-0.5*((Kh[0] - Kh[1]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[0] - Kh[1]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[0] - Kh[1]) . r0)])*
Exp[I*((Ks[0] - Ks[1] - Kh[0] + Kh[1]) . r)],
v1*(Exp[-0.5*((Kh[0] - Kh[2]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[0] - Kh[2]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[0] - Kh[2]) . r0)])*
Exp[I*((Ks[0] - Ks[2] - Kh[0] + Kh[2]) . r)],
v1*(Exp[-0.5*((Kh[0] - Kh[3]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[0] - Kh[3]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[0] - Kh[3]) . r0)])*
Exp[I*((Ks[0] - Ks[3] - Kh[0] + Kh[3]) . r)],
v1*(Exp[-0.5*((Kh[0] - Kh[4]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[0] - Kh[4]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[0] - Kh[4]) . r0)])*
Exp[I*((Ks[0] - Ks[4] - Kh[0] + Kh[4]) . r)],
v1*(Exp[-0.5*((Kh[0] - Kh[5]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[0] - Kh[5]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[0] - Kh[5]) . r0)])*
Exp[I*((Ks[0] - Ks[5] - Kh[0] + Kh[5]) . r)],
v1*(Exp[-0.5*((Kh[0] - Kh[6]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[0] - Kh[6]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[0] - Kh[6]) . r0)])*
Exp[I*((Ks[0] - Ks[6] - Kh[0] + Kh[6]) .
r)]}, {v1*(Exp[-0.5*((Kh[1] -
Kh[0]) . ((1 - Exp[(3 - rA*10)/0.3])*u))^2] +
Exp[-0.5*((Kh[1] - Kh[0]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[1] - Kh[0]) . r0)])*
Exp[I*((Ks[1] - Ks[0] - Kh[1] + Kh[0]) . r)], p[1],
v1*(Exp[-0.5*((Kh[1] - Kh[2]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[1] - Kh[2]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[1] - Kh[2]) . r0)])*
Exp[I*((Ks[1] - Ks[2] - Kh[1] + Kh[2]) . r)],
v2*(Exp[-0.5*((Kh[1] - Kh[3]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[1] - Kh[3]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[1] - Kh[3]) . r0)])*
Exp[I*((Ks[1] - Ks[3] - Kh[1] + Kh[3]) . r)],
v3*(Exp[-0.5*((Kh[1] - Kh[4]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[1] - Kh[4]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[1] - Kh[4]) . r0)])*
Exp[I*((Ks[1] - Ks[4] - Kh[1] + Kh[4]) . r)],
v2*(Exp[-0.5*((Kh[1] - Kh[5]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[1] - Kh[5]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[1] - Kh[5]) . r0)])*
Exp[I*((Ks[1] - Ks[5] - Kh[1] + Kh[5]) . r)],
v1*(Exp[-0.5*((Kh[1] - Kh[6]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[1] - Kh[6]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[1] - Kh[6]) . r0)])*
Exp[I*((Ks[1] - Ks[6] - Kh[1] + Kh[6]) .
r)]}, {v1*(Exp[-0.5*((Kh[2] -
Kh[0]) . ((1 - Exp[(3 - rA*10)/0.3])*u))^2] +
Exp[-0.5*((Kh[2] - Kh[0]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[2] - Kh[0]) . r0)])*
Exp[I*((Ks[2] - Ks[0] - Kh[2] + Kh[0]) . r)],
v1*(Exp[-0.5*((Kh[2] - Kh[1]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[2] - Kh[1]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[2] - Kh[1]) . r0)])*
Exp[I*((Ks[2] - Ks[1] - Kh[2] + Kh[1]) . r)], p[2],
v1*(Exp[-0.5*((Kh[2] - Kh[3]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[2] - Kh[3]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[2] - Kh[3]) . r0)])*
Exp[I*((Ks[2] - Ks[3] - Kh[2] + Kh[3]) . r)],
v2*(Exp[-0.5*((Kh[2] - Kh[4]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[2] - Kh[4]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[2] - Kh[4]) . r0)])*
Exp[I*((Ks[2] - Ks[4] - Kh[2] + Kh[4]) . r)],
v3*(Exp[-0.5*((Kh[2] - Kh[5]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[2] - Kh[5]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[2] - Kh[5]) . r0)])*
Exp[I*((Ks[2] - Ks[5] - Kh[2] + Kh[5]) . r)],
v2*(Exp[-0.5*((Kh[2] - Kh[6]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[2] - Kh[6]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[2] - Kh[6]) . r0)])*
Exp[I*((Ks[2] - Ks[6] - Kh[2] + Kh[6]) .
r)]}, {v1*(Exp[-0.5*((Kh[3] -
Kh[0]) . ((1 - Exp[(3 - rA*10)/0.3])*u))^2] +
Exp[-0.5*((Kh[3] - Kh[0]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[3] - Kh[0]) . r0)])*
Exp[I*((Ks[3] - Ks[0] - Kh[3] + Kh[0]) . r)],
v2*(Exp[-0.5*((Kh[3] - Kh[1]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[3] - Kh[1]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[3] - Kh[1]) . r0)])*
Exp[I*((Ks[3] - Ks[1] - Kh[3] + Kh[1]) . r)],
v1*(Exp[-0.5*((Kh[3] - Kh[2]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[3] - Kh[2]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[3] - Kh[2]) . r0)])*
Exp[I*((Ks[3] - Ks[2] - Kh[3] + Kh[2]) . r)], p[3],
v1*(Exp[-0.5*((Kh[3] - Kh[4]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[3] - Kh[4]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[3] - Kh[4]) . r0)])*
Exp[I*((Ks[3] - Ks[4] - Kh[3] + Kh[4]) . r)],
v2*(Exp[-0.5*((Kh[3] - Kh[5]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[3] - Kh[5]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[3] - Kh[5]) . r0)])*
Exp[I*((Ks[3] - Ks[5] - Kh[3] + Kh[5]) . r)],
v3*(Exp[-0.5*((Kh[3] - Kh[6]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[3] - Kh[6]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[3] - Kh[6]) . r0)])*
Exp[I*((Ks[3] - Ks[6] - Kh[3] + Kh[6]) .
r)]}, {v1*(Exp[-0.5*((Kh[4] -
Kh[0]) . ((1 - Exp[(3 - rA*10)/0.3])*u))^2] +
Exp[-0.5*((Kh[4] - Kh[0]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[4] - Kh[0]) . r0)])*
Exp[I*((Ks[4] - Ks[0] - Kh[4] + Kh[0]) . r)],
v3*(Exp[-0.5*((Kh[4] - Kh[1]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[4] - Kh[1]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[4] - Kh[1]) . r0)])*
Exp[I*((Ks[4] - Ks[1] - Kh[4] + Kh[1]) . r)],
v2*(Exp[-0.5*((Kh[4] - Kh[2]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[4] - Kh[2]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[4] - Kh[2]) . r0)])*
Exp[I*((Ks[4] - Ks[2] - Kh[4] + Kh[2]) . r)],
v1*(Exp[-0.5*((Kh[4] - Kh[3]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[4] - Kh[3]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[4] - Kh[3]) . r0)])*
Exp[I*((Ks[4] - Ks[3] - Kh[4] + Kh[3]) . r)], p[4],
v1*(Exp[-0.5*((Kh[4] - Kh[5]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[4] - Kh[5]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[4] - Kh[5]) . r0)])*
Exp[I*((Ks[4] - Ks[5] - Kh[4] + Kh[5]) . r)],
v2*(Exp[-0.5*((Kh[4] - Kh[6]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[4] - Kh[6]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[4] - Kh[6]) . r0)])*
Exp[I*((Ks[4] - Ks[6] - Kh[4] + Kh[6]) .
r)]}, {v1*(Exp[-0.5*((Kh[5] -
Kh[0]) . ((1 - Exp[(3 - rA*10)/0.3])*u))^2] +
Exp[-0.5*((Kh[5] - Kh[0]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[5] - Kh[0]) . r0)])*
Exp[I*((Ks[5] - Ks[0] - Kh[5] + Kh[0]) . r)],
v2*(Exp[-0.5*((Kh[5] - Kh[1]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[5] - Kh[1]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[5] - Kh[1]) . r0)])*
Exp[I*((Ks[5] - Ks[1] - Kh[5] + Kh[1]) . r)],
v3*(Exp[-0.5*((Kh[5] - Kh[2]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[5] - Kh[2]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[5] - Kh[2]) . r0)])*
Exp[I*((Ks[5] - Ks[2] - Kh[5] + Kh[2]) . r)],
v2*(Exp[-0.5*((Kh[5] - Kh[3]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[5] - Kh[3]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[5] - Kh[3]) . r0)])*
Exp[I*((Ks[5] - Ks[3] - Kh[5] + Kh[3]) . r)],
v1*(Exp[-0.5*((Kh[5] - Kh[4]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[5] - Kh[4]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[5] - Kh[4]) . r0)])*
Exp[I*((Ks[5] - Ks[4] - Kh[5] + Kh[4]) . r)], p[5],
v1*(Exp[-0.5*((Kh[5] - Kh[6]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[5] - Kh[6]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[5] - Kh[6]) . r0)])*
Exp[I*((Ks[5] - Ks[6] - Kh[5] + Kh[6]) .
r)]}, {v1*(Exp[-0.5*((Kh[6] -
Kh[0]) . ((1 - Exp[(3 - rA*10)/0.3])*u))^2] +
Exp[-0.5*((Kh[6] - Kh[0]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[6] - Kh[0]) . r0)])*
Exp[I*((Ks[6] - Ks[0] - Kh[6] + Kh[0]) . r)],
v1*(Exp[-0.5*((Kh[6] - Kh[1]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[6] - Kh[1]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[6] - Kh[1]) . r0)])*
Exp[I*((Ks[6] - Ks[1] - Kh[6] + Kh[1]) . r)],
v2*(Exp[-0.5*((Kh[6] - Kh[2]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[6] - Kh[2]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[6] - Kh[2]) . r0)])*
Exp[I*((Ks[6] - Ks[2] - Kh[6] + Kh[2]) . r)],
v3*(Exp[-0.5*((Kh[6] - Kh[3]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[6] - Kh[3]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[6] - Kh[3]) . r0)])*
Exp[I*((Ks[6] - Ks[3] - Kh[6] + Kh[3]) . r)],
v2*(Exp[-0.5*((Kh[6] - Kh[4]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[6] - Kh[4]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[6] - Kh[4]) . r0)])*
Exp[I*((Ks[6] - Ks[4] - Kh[6] + Kh[4]) . r)],
v1*(Exp[-0.5*((Kh[6] - Kh[5]) . ((1 - Exp[(3 - rA*10)/0.3])*
u))^2] +
Exp[-0.5*((Kh[6] - Kh[5]) . ((1 - Exp[(3 - rB*10)/0.3])*
u))^2]*Exp[-I*((Kh[6] - Kh[5]) . r0)])*
Exp[I*((Ks[6] - Ks[5] - Kh[6] + Kh[5]) . r)], p[6]}};
S1 = Transpose[Eigenvectors[M1]] .
DiagonalMatrix[
Exp[I*0.335*(Sqrt[(k*Cos[theta])^2 + Eigenvalues[M1]] -
k*Cos[theta])]] . Inverse[Transpose[Eigenvectors[M1]]];
S2 = Transpose[Eigenvectors[M2]] .
DiagonalMatrix[
Exp[I*0.335*(Sqrt[(k*Cos[theta])^2 + Eigenvalues[M2]] -
k*Cos[theta])]] . Inverse[Transpose[Eigenvectors[M2]]];
{{S11, S12, S13, S14, S15, S16, S17}, {S21, S22, S23, S24, S25,
S26, S27}, {S31, S32, S33, S34, S35, S36, S37}, {S41, S42, S43,
S44, S45, S46, S47}, {S51, S52, S53, S54, S55, S56, S57}, {S61,
S62, S63, S64, S65, S66, S67}, {S71, S72, S73, S74, S75, S76,
S77}} = S2 . S1;
T = S11*Conjugate[S11] + S21*Conjugate[S21] + S31*Conjugate[S31] +
S41*Conjugate[S41] + S51*Conjugate[S51] + S61*Conjugate[S61] +
S71*Conjugate[
S71] + (Conjugate[S11]*S21*Exp[I*((Kh[1] - Kh[0]) . r)] +
Conjugate[S11]*S31*Exp[I*((Kh[2] - Kh[0]) . r)] +
Conjugate[S11]*S41*Exp[I*((Kh[3] - Kh[0]) . r)] +
Conjugate[S11]*S51*Exp[I*((Kh[4] - Kh[0]) . r)] +
Conjugate[S11]*S61*Exp[I*((Kh[5] - Kh[0]) . r)] +
Conjugate[S11]*S71*Exp[I*((Kh[6] - Kh[0]) . r)] +
Conjugate[S21]*S31*Exp[I*((Kh[2] - Kh[1]) . r)] +
Conjugate[S21]*S41*Exp[I*((Kh[3] - Kh[1]) . r)] +
Conjugate[S21]*S51*Exp[I*((Kh[4] - Kh[1]) . r)] +
Conjugate[S21]*S61*Exp[I*((Kh[5] - Kh[1]) . r)] +
Conjugate[S21]*S71*Exp[I*((Kh[6] - Kh[1]) . r)] +
Conjugate[S31]*S41*Exp[I*((Kh[3] - Kh[2]) . r)] +
Conjugate[S31]*S51*Exp[I*((Kh[4] - Kh[2]) . r)] +
Conjugate[S31]*S61*Exp[I*((Kh[5] - Kh[2]) . r)] +
Conjugate[S31]*S71*Exp[I*((Kh[6] - Kh[2]) . r)] +
Conjugate[S41]*S51*Exp[I*((Kh[4] - Kh[3]) . r)] +
Conjugate[S41]*S61*Exp[I*((Kh[5] - Kh[3]) . r)] +
Conjugate[S41]*S71*Exp[I*((Kh[6] - Kh[3]) . r)] +
Conjugate[S51]*S61*Exp[I*((Kh[5] - Kh[4]) . r)] +
Conjugate[S51]*S71*Exp[I*((Kh[6] - Kh[4]) . r)] +
Conjugate[S61]*S71*Exp[I*((Kh[6] - Kh[5]) . r)]) +
Conjugate[
Conjugate[S11]*S21*Exp[I*((Kh[1] - Kh[0]) . r)] +
Conjugate[S11]*S31*Exp[I*((Kh[2] - Kh[0]) . r)] +
Conjugate[S11]*S41*Exp[I*((Kh[3] - Kh[0]) . r)] +
Conjugate[S11]*S51*Exp[I*((Kh[4] - Kh[0]) . r)] +
Conjugate[S11]*S61*Exp[I*((Kh[5] - Kh[0]) . r)] +
Conjugate[S11]*S71*Exp[I*((Kh[6] - Kh[0]) . r)] +
Conjugate[S21]*S31*Exp[I*((Kh[2] - Kh[1]) . r)] +
Conjugate[S21]*S41*Exp[I*((Kh[3] - Kh[1]) . r)] +
Conjugate[S21]*S51*Exp[I*((Kh[4] - Kh[1]) . r)] +
Conjugate[S21]*S61*Exp[I*((Kh[5] - Kh[1]) . r)] +
Conjugate[S21]*S71*Exp[I*((Kh[6] - Kh[1]) . r)] +
Conjugate[S31]*S41*Exp[I*((Kh[3] - Kh[2]) . r)] +
Conjugate[S31]*S51*Exp[I*((Kh[4] - Kh[2]) . r)] +
Conjugate[S31]*S61*Exp[I*((Kh[5] - Kh[2]) . r)] +
Conjugate[S31]*S71*Exp[I*((Kh[6] - Kh[2]) . r)] +
Conjugate[S41]*S51*Exp[I*((Kh[4] - Kh[3]) . r)] +
Conjugate[S41]*S61*Exp[I*((Kh[5] - Kh[3]) . r)] +
Conjugate[S41]*S71*Exp[I*((Kh[6] - Kh[3]) . r)] +
Conjugate[S51]*S61*Exp[I*((Kh[5] - Kh[4]) . r)] +
Conjugate[S51]*S71*Exp[I*((Kh[6] - Kh[4]) . r)] +
Conjugate[S61]*S71*Exp[I*((Kh[6] - Kh[5]) . r)]];
Re[k^2*Sin[2*theta]*T]];
G := NIntegrate[
integrand[x, y, theta, fi], {x, x0 - 0.5, x0 + 0.5}, {y, y0 - 0.5,
y0 + 0.5}, {theta, 0, Pi/2}, {fi, 0, 2*Pi}, AccuracyGoal -> 4,
PrecisionGoal -> 3, Method -> "MonteCarlo"];
G

revised version:
Clear["Global*"]
x0 = 0;
y0 = -50;
afa = Pi/180*0.2;
K0 = 29.5;
kc = 5.12317;
k = kc*Sqrt[5];
v1 = kc^2*4;
v2 = kc^2*2;
v3 = kc^2*1;
r0 = {0.142, 0, 0};
a1 = {0.213, 0.123, 0};
a2 = {0.213, -0.123, 0};
b1 = {14.75, 25.55, 0};
b2 = {14.75, -25.55, 0};
Kh[0] = {0, 0, 0};
Kh[1] = {K0, 0, 0};
Kh[2] = {K0/2, K0*Sqrt[3]/2, 0};
Kh[3] = {-K0/2, K0*Sqrt[3]/2, 0};
Kh[4] = {-K0, 0, 0};
Kh[5] = {-K0/2, -K0*Sqrt[3]/2, 0};
Kh[6] = {K0/2, -K0*Sqrt[3]/2, 0};
Khmat = DeveloperToPackedArray@
N[{{0., 0., 0.}, {K0, 0., 0.}, {K0/2, K0 Sqrt[3]/2, 0.}, {-K0/2,
K0 Sqrt[3]/2, 0.}, {-K0, 0., 0.}, {-K0/2, -K0 Sqrt[3]/2,
0.}, {K0/2, -K0 Sqrt[3]/2, 0.}}];
Ksmat = DeveloperToPackedArray@
N[{{0., 0., 0.}, {K0*Cos[afa], K0*Sin[afa],
0}, {K0*Cos[afa + Pi/3], K0*Sin[afa + Pi/3],
0}, {K0*Cos[afa + 2*Pi/3], K0*Sin[afa + 2*Pi/3],
0}, {K0*Cos[afa + Pi], K0*Sin[afa + Pi],
0}, {K0*Cos[afa + 4*Pi/3], K0*Sin[afa + 4*Pi/3],
0}, {K0*Cos[afa + 5*Pi/3], K0*Sin[afa + 5*Pi/3], 0}}];
V = ArrayFlatten[{{{{0.}},
ConstantArray[v1, {1, 6}]}, {ConstantArray[v1, {6, 1}],
ToeplitzMatrix[{0., v1, v2, v3, v2, v1}]}}];
M10 = V*(1. + Exp[I*Outer[Subtract, Khmat . r0, Khmat . r0, 1]]);
integrand[x_?NumericQ, y_, theta_, fi_] :=
Module[{ks, r, p, r1, x1, y1, r2, u, tranB, rB, tranA, rA, M20, M2,
eigs, U, S1, S2, S, T, S11, S21, S31, S41, S51, S61, S71},
ks = {k*Sin[theta]*Cos[fi], k*Sin[theta]*Sin[fi], k*Cos[theta]};
r = {x, y, 0};
p[i_] := k^2 - Total[(ks + Kh[i])^2];
r1 = {-y*afa, x*afa, 0};
x1 = Mod[r1 . b1/(2*Pi), 1];
y1 = Mod[r1 . b2/(2*Pi), 1];
r2 = x1*a1 + y1*a2;
u = First[SortBy[{r0 - r2, 2*r0 - r2}, Norm]];
tranB =
Min[Norm[r2], Norm[a1 - r2], Norm[a2 - r2], Norm[2*r0 - r2],
Norm[3*r0 - r2], Norm[a1 - r0 - r2], Norm[a2 - r0 - r2]];
rB = Sqrt[tranB^2 + 0.112];
tranA =
Min[Norm[r2], Norm[a1 - r2], Norm[a2 - r2], Norm[r0 - r2],
Norm[3*r0 - r2], Norm[a1 + r0 - r2], Norm[a2 + r0 - r2]];
rA = Sqrt[tranA^2 + 0.112];
pvec =
Table[k^2 - Dot[ks + Khmat[[i + 1]], ks + Khmat[[i + 1]]], {i, 1,
6}];
M1 = M10;
Do[M1[[i + 1, i + 1]] = p[i], {i, 1, 6}];
M20 = V*(Exp[-0.5*
Outer[Subtract, Khmat . ((1 - Exp[(3 - rA*10)/0.3])*u),
Khmat . ((1 - Exp[(3 - rA*10)/0.3])*u), 1]^2] +
Exp[-0.5*
Outer[Subtract, Khmat . ((1 - Exp[(3 - rB*10)/0.3])*u),
Khmat . ((1 - Exp[(3 - rB*10)/0.3])*u), 1]^2]*Exp[I*Outer[Subtract, Khmat . r0, Khmat . r0, 1]])*
Exp[I*Outer[Subtract, (Ksmat - Khmat) . r, (Ksmat - Khmat) . r,
1]];
M2 = M20;
Do[M2[[i + 1, i + 1]] = p[i], {i, 1, 6}];
{eigs, U} = Eigensystem[M1];
S1 = Transpose[U] .
DiagonalMatrix[
Exp[I*0.335*(Sqrt[(k*Cos[theta])^2 + eigs] - k*Cos[theta])]] .
Inverse[Transpose[U]];
{eigs1, U1} = Eigensystem[M2];
S2 = Transpose[U1] .
DiagonalMatrix[
Exp[I*0.335*(Sqrt[(k*Cos[theta])^2 + eigs1] - k*Cos[theta])]] .
Inverse[Transpose[U1]];
{S11, S21, S31, S41, S51, S61, S71} = (S2 . S1)[[All, 1]];
T = S11*Conjugate[S11] + S21*Conjugate[S21] + S31*Conjugate[S31] +
S41*Conjugate[S41] + S51*Conjugate[S51] + S61*Conjugate[S61] +
S71*Conjugate[
S71] + (Conjugate[S11]*S21*Exp[I*((Kh[1] - Kh[0]) . r)] +
Conjugate[S11]*S31*Exp[I*((Kh[2] - Kh[0]) . r)] +
Conjugate[S11]*S41*Exp[I*((Kh[3] - Kh[0]) . r)] +
Conjugate[S11]*S51*Exp[I*((Kh[4] - Kh[0]) . r)] +
Conjugate[S11]*S61*Exp[I*((Kh[5] - Kh[0]) . r)] +
Conjugate[S11]*S71*Exp[I*((Kh[6] - Kh[0]) . r)] +
Conjugate[S21]*S31*Exp[I*((Kh[2] - Kh[1]) . r)] +
Conjugate[S21]*S41*Exp[I*((Kh[3] - Kh[1]) . r)] +
Conjugate[S21]*S51*Exp[I*((Kh[4] - Kh[1]) . r)] +
Conjugate[S21]*S61*Exp[I*((Kh[5] - Kh[1]) . r)] +
Conjugate[S21]*S71*Exp[I*((Kh[6] - Kh[1]) . r)] +
Conjugate[S31]*S41*Exp[I*((Kh[3] - Kh[2]) . r)] +
Conjugate[S31]*S51*Exp[I*((Kh[4] - Kh[2]) . r)] +
Conjugate[S31]*S61*Exp[I*((Kh[5] - Kh[2]) . r)] +
Conjugate[S31]*S71*Exp[I*((Kh[6] - Kh[2]) . r)] +
Conjugate[S41]*S51*Exp[I*((Kh[4] - Kh[3]) . r)] +
Conjugate[S41]*S61*Exp[I*((Kh[5] - Kh[3]) . r)] +
Conjugate[S41]*S71*Exp[I*((Kh[6] - Kh[3]) . r)] +
Conjugate[S51]*S61*Exp[I*((Kh[5] - Kh[4]) . r)] +
Conjugate[S51]*S71*Exp[I*((Kh[6] - Kh[4]) . r)] +
Conjugate[S61]*S71*Exp[I*((Kh[6] - Kh[5]) . r)]) +
Conjugate[
Conjugate[S11]*S21*Exp[I*((Kh[1] - Kh[0]) . r)] +
Conjugate[S11]*S31*Exp[I*((Kh[2] - Kh[0]) . r)] +
Conjugate[S11]*S41*Exp[I*((Kh[3] - Kh[0]) . r)] +
Conjugate[S11]*S51*Exp[I*((Kh[4] - Kh[0]) . r)] +
Conjugate[S11]*S61*Exp[I*((Kh[5] - Kh[0]) . r)] +
Conjugate[S11]*S71*Exp[I*((Kh[6] - Kh[0]) . r)] +
Conjugate[S21]*S31*Exp[I*((Kh[2] - Kh[1]) . r)] +
Conjugate[S21]*S41*Exp[I*((Kh[3] - Kh[1]) . r)] +
Conjugate[S21]*S51*Exp[I*((Kh[4] - Kh[1]) . r)] +
Conjugate[S21]*S61*Exp[I*((Kh[5] - Kh[1]) . r)] +
Conjugate[S21]*S71*Exp[I*((Kh[6] - Kh[1]) . r)] +
Conjugate[S31]*S41*Exp[I*((Kh[3] - Kh[2]) . r)] +
Conjugate[S31]*S51*Exp[I*((Kh[4] - Kh[2]) . r)] +
Conjugate[S31]*S61*Exp[I*((Kh[5] - Kh[2]) . r)] +
Conjugate[S31]*S71*Exp[I*((Kh[6] - Kh[2]) . r)] +
Conjugate[S41]*S51*Exp[I*((Kh[4] - Kh[3]) . r)] +
Conjugate[S41]*S61*Exp[I*((Kh[5] - Kh[3]) . r)] +
Conjugate[S41]*S71*Exp[I*((Kh[6] - Kh[3]) . r)] +
Conjugate[S51]*S61*Exp[I*((Kh[5] - Kh[4]) . r)] +
Conjugate[S51]*S71*Exp[I*((Kh[6] - Kh[4]) . r)] +
Conjugate[S61]*S71*Exp[I*((Kh[6] - Kh[5]) . r)]];
Re[k^2*Sin[2*theta]*T]];
G := NIntegrate[
integrand[x, y, theta, fi], {x, x0 - 0.5, x0 + 0.5}, {y, y0 - 0.5,
y0 + 0.5}, {theta, 0, Pi/2}, {fi, 0, 2*Pi}, AccuracyGoal -> 4,
PrecisionGoal -> 3, Method -> "MonteCarlo"];
G

• Related questions (duplicates?): this and this. You seem to be asking the same question over and over. Commented Feb 11 at 18:35
• @Ghoster, No, they are not the same question. I have mentioned in this question:"A method has been previously provided by Ted Ersek. However, when there exists the code "u:= First[SortBy[{r0 - r2, 2*r0 - r2}, Norm]]", the method will give a wrong result because the function "First[SortBy[]]" is invalid." Commented Feb 12 at 0:38
• You reased the original question. Please keep that one at all time when the first answer has arrived. Otherwise things will get very confusing. You should add your edits at the bottom of your post, best with a descriptive title. These post are not meant as conversation, but as documents that are to archived, so that other people can search for them and leaarn from them. Commented Feb 12 at 13:16
• "(1) For Hint 1., I have builded M2 as your way, but the result is wrong (it can be verified by integrand[]). Can you help me to find the reason?" The problem here is that I have little idea what M2 is supposed to do. You made no attempt to explain it. And I simply do not have the time to reverse engineer what you want to do. My job here is to show some possibilities how you can make your code cleaner, less error prone, and faster rather than doing your programming work. So please do just what I would do in this situation: Try hard until it is correct. They call process this debugging. Commented Feb 12 at 13:25
• "(2) What is the meaning of "pvec=..."? It seems like can be deleted." Yes, that is indeed an artifact. I meant to replace $p[1]$, $p[2]$, ... by $pvec[[1]]$, $pvec[[2]]$, ... because indexing into arrays is a bit faster than a lookup in the downvalues of a symbol. But then I realized that there is only a handful of these calls, so it does not matter much. Commented Feb 12 at 13:27

There are many issues with the code. The main one is that it is looong and totally unmaintainable.

Hint 1.

I am pretty sure that that most of the matrix M1 is constant. So compute those parts only once and then add-in the parts that have to be update. Also, you really should write this as a loop or using some other constructs. Here a worked out example:

One-time code:

Khmat = DeveloperToPackedArray@N[{
{0., 0., 0.},
{K0, 0., 0.},
{K0/2, K0 Sqrt[3]/2, 0.},
{-K0/2, K0 Sqrt[3]/2, 0.},
{-K0, 0., 0.},
{-K0/2, -K0 Sqrt[3]/2, 0.},
{K0/2, -K0 Sqrt[3]/2, 0.}
}];

V = ArrayFlatten[{
{{{0.}}, ConstantArray[v1, {1, 6}]},
{ConstantArray[v1, {6, 1}],
ToeplitzMatrix[{0., v1, v2, v3, v2, v1}]}
}];

outer = Outer[Subtract, #, #, 1] &;

M10 = V (1. + Exp[I outer[Khmat.r0]]);


Now, in each call to integrand, evaluate

pvec = Table[ k^2 - Dot[ks + Khmat[[i + 1]], ks + Khmat[[i + 1]]], {i, 1, 6}];
M1 = M10;
Do[M1[[i + 1, i + 1]] = p[i], {i, 1, 6}];


This is way short and expresses your intent way better.

The matrix M2 can be build in a similar way.

Hint 2.

Use SetDelayed (:= only if you really want to delay something.

Quite likely, you want

r1 := {-y*afa, x*afa, 0};
x1 := Mod[r1 . b1/(2*Pi), 1];
y1 := Mod[r1 . b2/(2*Pi), 1];
r2 := x1*a1 + y1*a2;
u := First[SortBy[{r0 - r2, 2*r0 - r2}, Norm]];


but rather

r1 = {-y afa, x afa, 0};
x1 = Mod[r1 . b1/(2 Pi), 1];
y1 = Mod[r1 . b2/(2 Pi), 1];
r2 = x1 a1 + y1 a2;
u = First[SortBy[{r0 - r2, 2 r0 - r2}, Norm]];


Hint 3.

Use expensive computations like Eigenvalues and Eigenvectors only once, store their results. Moreover, you can fuse the calls to both of them by using Eigensystem. Moreover, Inverse[Transpose[Eigenvectors[M1]]] should be equivalent to Conjugate[Eigenvectors[M1]] for machine precision matrices due to the orthonormalization. So S1 can be computed as follows. (Mind that we can get rid of the DiagonalMatrix by using threaded multiplication.)

{eigs, U} = Eigensystem[M1];

S1 = Transpose[U].(Exp[I 0.335 (Sqrt[(k Cos[theta])^2 + eigs] - k Cos[theta])] Inverse[Transpose[Eigenvectors[M1]]]);


Hint 4.

The way you store the result of S1.S2 and use it afterwards just looks wrong. Also, you seem to use only the first column (S1.S2)[[All,1]] of the matrix-matrix product S1.S2, so you merely have to compute (S1.S2)[[All,1]] == S1.(S2[[All,1]]).

Edit

With regard to M20, you seem to have overlooked just one exponential factor. Also, you can use $$\mathrm{ee}^A \cdot \mathrm{ee}^A = \mathrm{ee}^{A+B}$$ to get rid of a couple of Exp evaluations. A call to arithmetic like + and * is several times faster than a call to Exp. I have not tested it (so no guarantees), but it appears to me that M20 should rather be constructed like this:

uA = ((1 - Exp[(3 - rA 10)/0.3]) u);
uB = ((1 - Exp[(3 - rB 10)/0.3]) u);

outer = Outer[Subtract, #, #, 1] &;

mat = I outer[(Ksmat - Khmat).r];

M20 = V (Exp[-0.5 outer[Khmat.uA]^2 + mat] +
Exp[-0.5 outer[Khmat.uB]^2 + I outer[Khmat.r0] + mat]);
`