% This code determines 2D bandgap structure in Fig. 7.5 for Eq. (7.39): 
% Uxx+Uyy+mu*U+V(x,y)U=0, where V(x, y)=-E0/(1+I0*cos(x)^2*cos(y)^2)

I0=6; E0=15; Lx=pi; Ly=pi; N=11;        % number of Fourier modes is 2N+1 
dx=Lx/100; dy=Ly/100; k0x=2*pi/Lx; k0y=2*pi/Ly; 
x=0:dx:Lx-dx; y=0:dy:Ly-dy; [X Y]=meshgrid(x, y); 
V=-E0./(1+I0*cos(X).^2.*cos(Y).^2);  
NK=20;     % NK is # of grid points along each segment of Gamma -> X -> M 
kxgrid=0:(pi/Lx)/NK:(pi/Lx); kygrid=0:(pi/Ly)/NK:(pi/Ly); 
Kpath=[kxgrid         pi/Lx*ones(1,NK)    kxgrid(NK:-1:1)  
       zeros(1,NK+1)  kygrid(2:NK+1)      kygrid(NK:-1:1)]; 

for m=1:1:2*N+1         % calculating Fourier coefficients of potential V
  for n=1:1:2*N+1
    [m n]
    integrand=V.*exp(-i*((m-N-1)*k0x*X+(n-N-1)*k0y*Y)); 
    Vfourier(m, n)=dx*dy*sum(sum(integrand))/(Lx*Ly); 
  end
end

for nn=1:length(Kpath)   % calculating Bloch curves along Gamma -> X -> M 
 nn
 k1=Kpath(1,nn); k2=Kpath(2,nn); L=zeros((2*N+1)^2); 
 for k=1:(2*N+1)^2
  for l=1:(2*N+1)^2
   ii=-N+floor((k-1)/(2*N+1));  j=-N+mod(k-1, 2*N+1);
    m=-N+floor((l-1)/(2*N+1));  n=-N+mod(l-1, 2*N+1);
   if ii-m+N+1 >= 1 & ii-m+N+1 <= 2*N+1 & j-n+N+1 >= 1 & j-n+N+1 <= 2*N+1 
      L(k, l)=-Vfourier(ii-m+N+1, j-n+N+1); 
   end   
   if k == l
      L(k, l)=L(k, l)+(k1+ii*k0x)^2+(k2+j*k0y)^2;   
   end
  end
 end
 mu=sort(eig(L)); Mu(nn,1:6)=mu(1:6);   
end

plot(0:3*NK, Mu(:,1:6)); hold on;                  % plotting Bloch bands 
plot([NK NK], [0 30], 'k'); plot([2*NK 2*NK], [0 30], 'k'); 
axis([0 3*NK 5 15]); ylabel('\mu', 'fontsize', 22); 
set(gca,'xtick', 0:NK:3*NK, 'xticklabel',{'','X','M',''}, 'fontsize',15); 
text(-NK/20, 4.5, '\Gamma', 'fontsize', 15)
text(3*NK-NK/20, 4.5, '\Gamma', 'fontsize', 15)