pw_te_tm.m

Calculating photonic bands with plane wave expansion method

clear all
warning off;
epsa=1.0;  %material for the atom
epsb=4.9912;    %material of the background
a=1.0;
R=0.3*a;  %radius of the atom
i=sqrt(-1);
f=pi*R^2/a^2;  %volume fraction of atom(s) in the cell
NumberofCell=1; %supercell size=1X1
a1=a;
a2=a*i;
b1=2*pi/a/NumberofCell;
b2=2*pi/a/NumberofCell*i;
n=5;%input('please input n:');
display('Fourier transforming.....');
tic;
NumberofPW=(2*n+1)^2;
mind=(-2*n:2*n)'+2*n+1;
mind=mind(:,ones(1,size(mind)))-2*n-1;
GG=mind'*b1+mind*b2;
%clear mind;
eps21=2*f*(epsa-epsb)*besselj(1,abs(GG).*R)./(abs(GG).*R);
eps21(2*n+1,2*n+1)=epsb+f*(epsa-epsb);
zz=[0,0]*[a1 a2].';%use 1X1 supercell to verify the algorithm
eps20=zeros(length(eps21));
for x=1:length(zz),
   eps20=eps20+exp(i*(real(GG).*real(zz(x))+imag(GG).*imag(zz(x)))).*eps21;
end
ff=pi*R^2*length(zz)/(NumberofCell^2*a^2);
eps20=eps20./NumberofCell^2;
eps20(2*n+1,2*n+1)=epsb+ff*(epsa-epsb);
count=1;
for y=-n:n,
 for x=-n:n;
      G(count)=x*b1+y*b2;
      r(count,:)=[x,y];
      count=count+1;
   end
end

display('Building eps(G,G) matrix from the FFT matrix');
for x=1:NumberofPW,
  	for y=x:NumberofPW,
	b=r(x,:)-r(y,:)+2*n+1;        
	eps2(x,y)=eps20(b(2),b(1));
        eps2(y,x)=eps2(x,y);   
   end
end
   

k2=2*pi/a*0.5.*(0.1:0.1:1.0);           %Gamma-X
k3=2*pi/a*(0.5+(0.0:0.1:1.0).*0.5*i);   %X-M
%c1=(1-0)/(sqrt(2)*length(k2));
k1=2*pi/a*(0.5+0.5*i).*(1.0:-0.1:0.0);     %M-Gamma

k0=[k1 k2 k3];%./NumberofCell;

% k0 = k1;

counter=1;
eps2=inv(eps2);

for ii=1:length(k0),
    k=k0(ii);
    % M=(real(k+G.')*real(k+G)+imag(k+G.')*imag(k+G)).*eps2; %TE wave
    M=abs(k+G.')*abs(k+G).*eps2; %TM wave
	E=sort(abs(eig(M))); 
    freq(:,counter)=sqrt(abs(E)).*a./2./pi; 
	display(sprintf('calculation of k=%f+%fi is finished',real(k),imag(k)));
	counter=counter+1;
end
toc;
tmpx=1:length(k0);
plot(tmpx,freq,'k','linewidth',1); 
%title('E-polarization, f=0.3, r/a=0.3090');
xlabel('wave vector')
ylabel('Normalized frequency \omegaa/2\pic')
grid on

axis([1 length(k0) 0 0.6])
y1=get(gca,'ylim');
x1=get(gca, 'xlim');
kl1=[length(k1) length(k1)];
kl2=[length(k1)+length(k2) length(k1)+ length(k2)];

line(kl1,y1,'Color','k');
line(kl2,y1,'Color','k');

b=min(ylim);
c=b-0.04;
h1=text(min(xlim),c,'M');
h2=text(length(k1),c,'\Gamma');
h3=text(length(k1)+length(k2),c,'X');
h4=text(length(k0),c,'M');
% h3=text(length(k3)+length(k2)+length(k1),c,'M');
ax=gca;


ax=gca;
set(gca,'XTick',[]);
% set(ax,'XTickLabel',[0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5]);
set(ax,'XTickLabel',[]);
set(gca,'YTick',[0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 1.05 1.2]);
set(gca,'YTickLabel',[0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 1.05 1.2])