goalp=2
period=2 π
kmax=(5/2)
kperpmax = kmax
kmin=(1/128)
nmax[X]=Floor[((kperpmax period)/(2 π))]
nmin[X]=Ceiling[((-kperpmax period)/(2 π))]
nmax[Y]=0
nmin[Y]=0
kperp[kx0_,nx_,X]:=((2 π nx)/period)+kx0
kperp[ky0_,0,Y]:=ky0
kperpmagsq[kx0_,ky0_,nx_,ny_]:=((kperp[kx0,nx,X])^2)+((kperp[ky0,ny,Y])^2)
kz[k_,kx0_,ky0_,nx_,ny_]:=k(Sqrt[1-((kperpmagsq[kx0,ky0,nx,ny])/(k^2))])
allowedWave[nx_,ny_]:=nx≤nmax[X]&& nx≥nmin[X] && ny ≤nmax[Y]&&ny ≥nmin[Y]
indexLookupTable=Module[{nx,ny,table={}},For[nx=nmin[X],nx≤nmax[X],nx++,
For[ny=nmin[Y],ny≤nmax[Y],ny++,If[allowedWave[nx,ny],table=Append[table,{nx,ny}]]]];table]

maxIndex=3Length[indexLookupTable]
indexton[i_,comp_]:=indexLookupTable[[Floor[(i-1)/3]+1]][[If[comp===X,1,2]]]
indextocomp[i_]:=Switch[Mod[i-1,3],0,X,1,Y,2,Z]

DOperator[expr_]:=D[expr,z]

gradOperator[(A_:1)(L1_+L2_)]:=gradOperator[A L1]+gradOperator[A L2]

gradOperator[(L1_:1) planeWave[k_,kx0_,ky0_,nx_,ny_]]/;(FreeQ[L1,planeWaveVector]&&FreeQ[L1,planeWave]):=I kperp[kx0,nx,X]L1 planeWaveVector[k,kx0,ky0,nx,ny,X]+I kperp[ky0,ny,Y]L1 planeWaveVector[k,kx0,ky0,nx,ny,Y]+DOperator[L1] planeWaveVector[k,kx0,ky0,nx,ny,Z]

gradOperator[expr_]/;(FreeQ[expr,planeWave]&&FreeQ[expr,planeWaveVector]&&FreeQ[expr,z]):=0

divOperator[(A_:1)(L1_+L2_)]:=divOperator[A L1]+divOperator[A L2]

divOperator[(L1_:1)planeWaveVector[k_,kx0_,ky0_,nx_,ny_,comp_]]/;(FreeQ[L1,planeWaveVector]&&FreeQ[L1,planeWave]):=
Switch[comp,X,I kperp[kx0,nx,X]L1 planeWave[k,kx0,ky0,nx,ny],Y,I kperp[ky0,ny,Y]L1 planeWave[k,kx0,ky0,nx,ny],Z,DOperator[L1] planeWave[k,kx0,ky0,nx,ny]]

divOperator[L1_]/;(FreeQ[L1,z]&&FreeQ[expr,planeWave]&&FreeQ[L1,planeWaveVector]):=0

curlOperator[(A_:1)(L1_+L2_)]:=curlOperator[A L1]+curlOperator[A L2]

curlOperator[(L1_:1)planeWaveVector[k_,kx0_,ky0_,nx_,ny_,comp_]]/;(FreeQ[L1,planeWaveVector]&&FreeQ[L1,planeWave]):=
Switch[comp,X,DOperator[L1]planeWaveVector[k,kx0,ky0,nx,ny,Y]-I kperp[ky0,ny,Y]L1 planeWaveVector[k,kx0,ky0,nx,ny,Z],Y,-DOperator[L1]planeWaveVector[k,kx0,ky0,nx,ny,X]+I kperp[kx0,nx,X]L1 planeWaveVector[k,kx0,ky0,nx,ny,Z],Z,I kperp[ky0,ny,Y]L1 planeWaveVector[k,kx0,ky0,nx,ny,X]-I kperp[kx0,nx,X]L1 planeWaveVector[k,kx0,ky0,nx,ny,Y]]

(*This piece is necessary mostly so that curlOperator[0]=0. Any other argument should be some constant or function of r times a vector plane wave.*)
curlOperator[L1_]/;(FreeQ[L1,z]&&FreeQ[expr,planeWave]&&FreeQ[L1,planeWaveVector]):=0

(*This operator is equivalent to multiplying (dot product) by the complex conjugate of planeWaveVector[k,kx0,n,comp] and integrating over the perpendicular directions *)

perpIntegrate[(A_:1)(L1_+L2_),nx_,ny_,comp_]:=perpIntegrate[A L1,nx,ny,comp]+perpIntegrate[A L2,nx,ny,comp]

(*If the angular dependent part is not a PWV, perpIntegrate gives zero. This is not a problem if everything is expanded as PWV*)
perpIntegrate[A_ B_,nx_,ny_,comp_]/;(FreeQ[A,planeWaveVector]&&FreeQ[A,planeWave]):=A perpIntegrate[B,nx,ny,comp]
perpIntegrate[planeWaveVector[k_,kx0_,ky0_,nx1_,ny1_,comp1_],nx2_,ny2_,comp2_]:=If[comp1===comp2,KroneckerDelta[nx1,nx2]KroneckerDelta[ny1,ny2],0]

(*A "product" function to multiply scalar and vector waves together.*)

product[planeWave[k_,kx0_,ky0_,nx1_,ny1_],planeWaveVector[k_,kx0_,ky0_,nx2_,ny2_,comp_]]:=If[allowedWave[nx1+nx2,ny1+ny2],planeWaveVector[k,kx0,ky0,nx1+nx2,ny1+ny2,comp],0]

product[L1_+L2_,L3_]:=product[L1,L3]+product[L2,L3]
product[L1_,L2_+L3_]:=product[L1,L2]+product[L1,L3]
product[A_ L1_,L2_]/;(FreeQ[A,planeWave]):=A product[L1,L2]
product[L1_,A_ L2_]/;(FreeQ[A,planeWaveVector]):=A product[L1,L2]

evector[k_,kx0_,ky0_,z_]:=Sum[(Subscript[e, {indexton[i,X],indexton[i,Y],indextocomp[i]}]) [z]planeWaveVector[k,kx0,ky0,indexton[i,X],indexton[i,Y],indextocomp[i]],{i,1,maxIndex}]

myδ[nx_,ny_,comp_,ip_]:=If[comp===indextocomp[ip],KroneckerDelta[nx,indexton[ip,X]]KroneckerDelta[ny,indexton[ip,Y]],0]

getkcomp[k_,kx0_,ky0_,i_]:=Switch[indextocomp[i],X, kperp[kx0,indexton[i,X],X],Y,kperp[ky0,indexton[i,Y],Y],Z,kz[k,kx0,ky0,indexton[i,X],indexton[i,Y]]]

width=2
steepness=8
depth[k_]=2
profile[k_,z_]:=depth[k](1/(1+Exp[2 steepness (z-width)]))

ϵ[k_,nx_,ny_]/;(ny≠0||Abs[nx]>1):=Function[0]

ϵtotal[k_,kx0_,ky0_,z_]:=Sum[ϵ[k,nx,ny][z] planeWave[k,kx0,ky0,nx,ny],{nx,nmin[X],nmax[X]},{ny,nmin[Y],nmax[Y]}]

DistributeDefinitions["Global`"]
ParallelTable[{$RecursionLimit=∞,$MaxExtraPrecision=1024},{i,1,$KernelCount}]
$RecursionLimit=∞;
$MaxExtraPrecision=1024
SessionTime[]
system[k_,kx0_,ky0_]=Module[{curlcurl,ϵevec,ϵdivϵ,eqn},
{curlcurl,ϵevec}=Parallelize[{curlOperator[curlOperator[evector[k,kx0,ky0,z]]],
product[ϵtotal[k,kx0,ky0,z],evector[k,kx0,ky0,z]]}];
ϵdivϵ=product[ϵtotal[k,kx0,ky0,z],gradOperator[divOperator[ϵevec]]];
eqn=Expand[curlcurl-ϵdivϵ-(k^2)ϵevec];
DistributeDefinitions[eqn];
ParallelTable[perpIntegrate[eqn,indexton[i,X],indexton[i,Y],indextocomp[i]],{i,1,maxIndex}]
];
SessionTime[]
zerorule={x_/;x==0->0}
getd02d12[k_,kx0_,ky0_]=Module[{sys=system[k,kx0,ky0],d,d1,d2,d2inv},
DistributeDefinitions[sys,myδ,maxIndex];{d2,d1,d}=Parallelize[{
Table[FullSimplify[-sys[[i]]/.{(Subscript[e, {nx_,ny_,comp_}])[z]->0,(((Subscript[e, {nx_,ny_,comp_}])))'[z]->0,(((Subscript[e, {nx_,ny_,comp_}])))''[z]:>myδ[nx,ny,comp,ip]}],{i,1,maxIndex},{ip,1,maxIndex}],
Table[FullSimplify[sys[[i]]/.{(Subscript[e, {nx_,ny_,comp_}])[z]->0,(((Subscript[e, {nx_,ny_,comp_}])))''[z]->0,(((Subscript[e, {nx_,ny_,comp_}])))'[z]:>myδ[nx,ny,comp,ip]}],{i,1,maxIndex},{ip,1,maxIndex}],
Table[FullSimplify[sys[[i]]/.{(((Subscript[e, {nx_,ny_,comp_}])))'[z]->0,(((Subscript[e, {nx_,ny_,comp_}])))''[z]->0,(Subscript[e, {nx_,ny_,comp_}])[z]:>myδ[nx,ny,comp,ip]}],{i,1,maxIndex},{ip,1,maxIndex}]}];
d2inv=FullSimplify[Inverse[d2,Method->OneStepRowReduction]];
DistributeDefinitions[d2inv,d,d1];
Parallelize[{FullSimplify[d2inv.d],FullSimplify[d2inv.d1]}]/.zerorule];
{getd02A[k_,kx0_,ky0_],getd12A[k_,kx0_,ky0_]}=getd02d12[k,kx0,ky0];
Clear[getd02d12]

getd12PrimeA[k_,kx0_,ky0_]=FullSimplify[D[getd12A[k,kx0,ky0],z]];
SessionTime[]
ϵ[k_,0,0]:=Function[z,profile[k,z]+1]
ϵ[k_,1,0]:=Function[z,(1/2)profile[k,z]]
ϵ[k_,-1,0]:=Function[z,(1/2)profile[k,z]]
DistributeDefinitions[ϵ,profile,getd02A, getd12A, getd12PrimeA];
{getd02[k_,kx0_,ky0_],getd12[k_,kx0_,ky0_],getd12Prime[k_,kx0_,ky0_]}=Parallelize[{Simplify[Factor[Simplify[getd02A[k,kx0,ky0]]]],Simplify[Factor[Simplify[getd12A[k,kx0,ky0]]]],Simplify[Factor[Simplify[getd12PrimeA[k,kx0,ky0]]]]}];
getd2diff[k_,kx0_,ky0_]=Simplify[getd02[k,kx0,ky0]-getd12Prime[k,kx0,ky0]];
Clear[getd02A, getd12A, getd12PrimeA];
ParallelTable[Clear[getd02A, getd12A, getd12PrimeA],{i,1,$KernelCount}]
SessionTime[]
free[k_,kx0_,ky0_,nx_,ny_,z_]:=Exp[I kz[k,kx0,ky0,nx,ny]z]
logderivFree[k_,kx0_,ky0_,nx_,ny_]:=I kz[k,kx0,ky0,nx,ny]
compSign[i_,anti]:=If[indextocomp[i]===Z,-1,1]
compSign[i_,sym]:=If[indextocomp[i]===Z,1,-1]

freeMatrix[k_,kx0_,ky0_,z_]:=DiagonalMatrix[Table[free[k,kx0,ky0,indexton[i,X],indexton[i,Y],z],{i,1,maxIndex}]]
freeMatrixInverseSign[k_,kx0_,ky0_,z_,symanti_]:=DiagonalMatrix[Table[(1/(free[compSign[i,symanti]k,kx0,ky0,indexton[i,X],indexton[i,Y],z])),{i,1,maxIndex}]]

logderivFreeMatrix[k_,kx0_,ky0_]:=DiagonalMatrix[Table[logderivFree[k,kx0,ky0,indexton[i,X],indexton[i,Y]],{i,1,maxIndex}]]
fluxnorm[k_,kx0_,ky0_]:=DiagonalMatrix[Table[(Sqrt[(k/(kz[k,kx0,ky0,indexton[i,X],indexton[i,Y]]))]),{i,1,maxIndex}]]
kzsqMatrix[k_,kx0_,ky0_]:=DiagonalMatrix[Table[((kz[k, kx0,ky0,indexton[i,X],indexton[i,Y]])^2),{i,1,maxIndex}]]

signMatrix[symanti_] := DiagonalMatrix[Table[compSign[i,symanti],{i,1,maxIndex}]]
hzeroMatrix[k_,kx0_,ky0_,symanti_]:=DiagonalMatrix[Table[Switch[compSign[i,symanti],-1,(1/(logderivFree[-k,kx0,ky0,indexton[i,X],indexton[i,Y]])),1,0],{i,1,maxIndex}]]
longitudinalVector[k_,kx0_,ky0_,nx_,ny_]:=Table[KroneckerDelta[indexton[i,X],nx]KroneckerDelta[indexton[i,Y],ny]((getkcomp[k,kx0,ky0,i])/k),{i,1,maxIndex}]

transverseVectorTE[k_,kx0_,ky0_,nx_,ny_]:=Table[KroneckerDelta[indexton[i,X],nx]KroneckerDelta[indexton[i,Y],ny]If[kperpmagsq[kx0,ky0,nx,ny]==0, If[indextocomp[i]===Y,1,0],
((Switch[indextocomp[i],X,kperp[ky0,ny,Y],Y,-kperp[kx0,nx,X],Z,0])/(Sqrt[kperpmagsq[kx0,ky0,nx,ny]]))],{i,1,maxIndex}]
transverseVectorTM[k_,kx0_,ky0_,nx_,ny_]:=Table[KroneckerDelta[indexton[i,X],nx]KroneckerDelta[indexton[i,Y],ny]If[kperpmagsq[kx0,ky0,nx,ny]==0, If[indextocomp[i]===X,1,0],
((getkcomp[k,kx0,ky0,i]kz[k,kx0,ky0,nx,ny]-If[indextocomp[i]===Z,(k^2),0])/(k(Sqrt[kperpmagsq[kx0,ky0,nx,ny]])))],{i,1,maxIndex}]

infinity=2 width
zero=0
matrixsystemg[k_,kx0_,ky0_]:=-gmatrix''[z]+(getd12[k,kx0,ky0].gmatrix[z]-2 gmatrix'[z]).logderivFreeMatrix[k,kx0,ky0]+getd12[k,kx0,ky0].gmatrix'[z]+getd02[k,kx0,ky0].gmatrix[z]+gmatrix[z].kzsqMatrix[k,kx0,ky0]==0
boundary1g[k_,kx0_,ky0_]:=gmatrix[infinity]==IdentityMatrix[maxIndex]
boundary2g[k_,kx0_,ky0_]:=gmatrix'[infinity]==DiagonalMatrix[Table[0,{i,1,maxIndex}]]

matrixsystemh[k_,kx0_,ky0_,symanti_]:=-hmatrix''[z]+signMatrix[symanti].logderivFreeMatrix[k,kx0,ky0].(2 hmatrix'[z]+hmatrix[z].getd12[k,kx0,ky0])-hmatrix'[z].getd12[k,kx0,ky0]+hmatrix[z].getd2diff[k,kx0,ky0]+kzsqMatrix[k,kx0,ky0].hmatrix[z]==0
boundary1h[k_,kx0_,ky0_,symanti_]:=hmatrix'[zero]==IdentityMatrix[maxIndex]
boundary2h[k_,kx0_,ky0_,symanti_]:=hmatrix[zero]==hzeroMatrix[k,kx0,ky0,symanti]

solg[k_,kx0_,ky0_,x_:zero]:=First[NDSolve[{matrixsystemg[k,kx0,ky0],boundary1g[k,kx0,ky0],boundary2g[k,kx0,ky0]},gmatrix,{z,x,infinity},MaxSteps->100000,Method->StiffnessSwitching]]
solh[k_,kx0_,ky0_,x_:infinity,symanti_]:=First[NDSolve[{matrixsystemh[k,kx0,ky0,symanti],boundary1h[k,kx0,ky0,symanti],boundary2h[k,kx0,ky0,symanti]},hmatrix,{z,zero,x},MaxSteps->100000,Method->StiffnessSwitching]]

wronk[g_,h_,x_,k_,kx0_,ky0_,symanti_]:=freeMatrixInverseSign[k,kx0,ky0,x,symanti].(h[x].(g[x].logderivFreeMatrix[k,kx0,ky0]+g'[x])-(h'[x]-signMatrix[symanti].logderivFreeMatrix[k,kx0,ky0].h[x]).g[x] - h[x].(getd12[k,kx0,ky0]/.{z->x}).g[x]).freeMatrix[k,kx0,ky0,x].fluxnorm[k,kx0,ky0]

smatdiff[k_,kx0_,ky0_,x_]:=Module[{g,gm,hs,hms,ha,hma},
{g,gm, hs,hms,ha,hma}=Parallelize[{
gmatrix/.solg[k,kx0,ky0,x],gmatrix/.solg[-k,kx0,ky0,x],hmatrix/.solh[k,kx0,ky0,x,sym],hmatrix/.solh[-k,kx0,ky0,x,sym],hmatrix/.solh[k,kx0,ky0,x,anti],hmatrix/.solh[-k,kx0,ky0,x,anti]}];(Inverse[wronk[g,hs,x,k,kx0,ky0,sym]].signMatrix[sym].wronk[gm,hms,x,-k,kx0,ky0,sym].signMatrix[sym]-Inverse[wronk[g,ha,x,k,kx0,ky0,anti]].signMatrix[anti].wronk[gm,hma,x,-k,kx0,ky0,anti].signMatrix[anti])]

makeVectorList[k_,kx0_,ky0_,type_]:=Module[{nx,ny,table={}},For[nx=nmin[X],nx≤nmax[X],nx++,
For[ny=nmin[Y],ny≤nmax[Y],ny++,If[allowedWave[nx,ny],table=Append[table,type[k,kx0,ky0,nx,ny]]]]];table]

gettetmlonlist[k_,kx0_,ky0_]:={makeVectorList[k,kx0,ky0,transverseVectorTE],makeVectorList[k,kx0,ky0,transverseVectorTM],makeVectorList[k,kx0,ky0,longitudinalVector]}

projectMatrix[a_,k_,kx0_,ky0_]:=Module[{u,telist, tmlist, lonlist,picklist},
{telist,tmlist,lonlist}=gettetmlonlist[k,kx0,ky0];
u=Join[telist,tmlist,lonlist];
picklist=Join[Table[True,{Length[telist]+Length[tmlist]}],Table[False,{Length[lonlist]}]];
Pick[Map[Pick[#,picklist]&,u.a.Transpose[u]],picklist]]

rmatp[k_,kx0_,ky0_,x_]:=(1/2)projectMatrix[smatdiff[k,kx0,ky0,x],k,kx0,ky0]

translation[k_,kx0_,ky0_,sep_,shift_]:=projectMatrix[DiagonalMatrix[Table[Exp[I kz[k,kx0,ky0,indexton[i,X],indexton[i,Y]]sep]Exp[I kperp[kx0,indexton[i,X],X]shift],{i,1,maxIndex}]],k,kx0,ky0]

integrand[k_?NumberQ,kx0_?NumberQ,ky0_?NumberQ,sep_,shift_]:=Module[{r=rmatp[k,kx0,ky0,(width/2)],tr1=translation[k,kx0,ky0,sep,shift],tr2=translation[k,kx0,ky0,sep,-shift],val},val=Log[Det[IdentityMatrix[((2maxIndex)/3)]-tr1.r.tr2.r]];
Print[k," ",kx0," ",ky0," ",val," ",SessionTime[]]; val]

DistributeDefinitions["Global`"]

sep=6
sh=(period/4)

SessionTime[]
AbsoluteTiming[NIntegrate[(1/(2(π^3)))integrand[I κ,kx,ky,sep,sh],{κ,kmin,kmax},{kx,0,(π/period)},{ky,kmin,kmax},PrecisionGoal->goalp]]

NotebookSave[]