findmode_kz.m

两个平面波展开程序

% this program is designed for solving the two dimensional photonic structure
% Photonic Crystal Fiber: Band gap with a kz
% Triangular lattice
% Shangping Guo
% Old Dominion University, 2000
clear
warning off
ea=1;
eb=1.46^2;  % The dielectric constant of background
R = 0.48;
PCType = 2;
kz_const = 1.0;
%define the lattice vectors
a=1.0;  % lattice constance
if PCType == 1  %Square lattice
    a1 = a*[1,0,0];
    a2 = a*[0,1,0];
    a3 = a*[0,0,1];
end
if PCType == 2  % Triangular lattice
    a1 = a*[1/2,sqrt(3)/2,0];
    a2 = a*[1/2,-sqrt(3)/2,0];
    a3 = a*[0,0,1];
end

%calculate the reciprocal lattice vectors
ac=abs(a1(1)*a2(2)-a1(2)*a2(1)); 
f = pi*R^2/ac;
volcell=dot(a1,cross(a2,a3));
b1 = 2*pi*cross(a2,a3)/volcell;
b2 = 2*pi*cross(a3,a1)/volcell;
b3 = 2*pi*cross(a1,a2)/volcell;

MaxDimForG = 10;    % The max positive number of the reciprocal lattice, G
DimForG = 2*MaxDimForG + 1;
NPW=DimForG^2;  % NPW: The number of Plane Waves

%
%            ^ 
%            | Y
%      O  O  O  O  O  -->(MaxDimForG,MaxDimForG)for this point!
%      O  O  O  O  O 
%    --O--O--O--O--O--> X
%      O  O  O  O  O
%      O  O  O  O  O 
%            |   
%            |

gtemp=-MaxDimForG:MaxDimForG;
gtemp1=repmat(gtemp,DimForG,1);
Gx=b1(1)*gtemp1+b2(1)*gtemp1';
Gy=b1(2)*gtemp1+b2(2)*gtemp1';
Gx=Gx(:)';
Gy=Gy(:)';

disp(strcat('The number of plane waves is--',num2str(NPW)));

Gx_m=repmat(Gx,NPW,1);
Gx_n=Gx_m';
Gy_m=repmat(Gy,NPW,1);
Gy_n=Gy_m'; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% count=1;
% for x = -MaxDimForG:MaxDimForG
%   	for y = -MaxDimForG:MaxDimForG
%         G(count,:) = x*b1 + y*b2;
%         count = count + 1;
%     end
% end
% 
% eps2 = zeros(NumberofPW,NumberofPW);
% for x=1:NumberofPW,
%   	for y=x+1:NumberofPW,
%    	eps2(x,y)=(epsa-epsb)*2*f*besselj(1,norm(G(x,:)-G(y,:))*R)./(norm(G(x,:)-G(y,:))*R);
%    	eps2(y,x)=eps2(x,y);   
%    end
%   	eps2(x,x)=f*epsa+(1-f)*epsb;
% end
% eps2=inv(eps2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 1024;
ic = N/2;
jc = N/2;
eps1 = zeros(N,N);
xdist_temp = zeros(1,5);
ydist_temp = zeros(1,5);
dist_temp = zeros(1,5);
for j = 1:N
    for i = 1:N
        xdist_temp(1) = (i - ic)*a1(1) + (j - jc)*a2(1);
        ydist_temp(1) = (i - ic)*a1(2) + (j - jc)*a2(2);
        dist_temp(1) = sqrt((xdist_temp(1)/N)^2 + (ydist_temp(1)/N)^2);
        xdist_temp(2) = (i - ic - N)*a1(1) + (j - jc)*a2(1);
        ydist_temp(2) = (i - ic - N)*a1(2) + (j - jc)*a2(2);
        dist_temp(2) = sqrt((xdist_temp(2)/N)^2 + (ydist_temp(2)/N)^2);
        xdist_temp(3) = (i - ic + N)*a1(1) + (j - jc)*a2(1);
        ydist_temp(3) = (i - ic + N)*a1(2) + (j - jc)*a2(2);
        dist_temp(3) = sqrt((xdist_temp(3)/N)^2 + (ydist_temp(3)/N)^2);
        xdist_temp(4) = (i - ic)*a1(1) + (j - jc - N)*a2(1);
        ydist_temp(4) = (i - ic)*a1(2) + (j - jc - N)*a2(2);
        dist_temp(4) = sqrt((xdist_temp(4)/N)^2 + (ydist_temp(4)/N)^2);
        xdist_temp(5) = (i - ic)*a1(1) + (j - jc + N)*a2(1);
        ydist_temp(5) = (i - ic)*a1(2) + (j - jc + N)*a2(2);
        dist_temp(5) = sqrt((xdist_temp(5)/N)^2 + (ydist_temp(5)/N)^2);
        dist = min(dist_temp);
        if (dist <= R)
            eps1(i,j) = 1/ea;
        else
            eps1(i,j) = 1/eb;
        end
    end
end
eps1_shift = fftshift(eps1);
eps20 = fftshift(fft2(eps1_shift)/N^2);
ek_mat = zeros(NPW,NPW);
for jj = 1:NPW
    ic_new = ic + floor((jj-1)/DimForG)+1;
    jc_new = jc + rem((jj-1),DimForG)+1;
    for ii = 1:NPW
        ii_new = ic_new - floor((ii-1)/DimForG);
        jj_new = jc_new - rem((ii-1),DimForG);
        ek_mat(ii,jj) = eps20(ii_new,jj_new);
    end
end
ek_mat = real(ek_mat);
%figure,mesh(ek_mat);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %k1=(0:0.1:1.0)/sqrt(3).*i*2*pi/a;
% %k2=((0.1:0.1:1.0)./3+1/sqrt(3)*i).*2.*pi./a;
% %k3=(0.9:-0.1:0).*(1.0/3.0+1/sqrt(3)*i).*2.*pi./a;%-(1/3+1/sqrt(3)*i)*2*pi/a;
% %k0=[k1 k2 k3];
% %k0=k1;
% %k0=(1/3+1/sqrt(3)*i).*2*pi/a;
% k0=zeros(length(0:0.1:1)+length(0.1:0.1:1)+length(0.9:-0.1:0),3);
% mm0=length(0:0.1:1);
% k0(1:mm0,2)=2*pi/a/sqrt(3).*(0:0.1:1)';
% mm0=length(0:0.1:1);
% mm1=length(0.1:0.1:1);
% k0(mm0+1:mm0+mm1,1)=2*pi/a/3.*(0.1:0.1:1)';
% k0(mm0+1:mm0+mm1,2)=2*pi/a/sqrt(3);
% mm0=mm0+mm1+1;
% mm1=length(0.9:-0.1:0);
% k0(mm0:mm0+mm1-1,1)=2*pi/a/3.*(0.9:-0.1:0)';
% k0(mm0:mm0+mm1-1,2)=2*pi/a/sqrt(3).*(0.9:-0.1:0)';
% k0(:,3)=9.0;
% 
% %k0=[1/3 1/sqrt(3) 9]*2*pi/a;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ktype = 3;  % The number of band parts, such as Gamma->K, K->M,...
Keach = 10; % The number of k vectors in each wave vector branch.
NumberK = Ktype*Keach;  % The total number of K vector.
NEIG = 20;  % The cut off number to the eigenvalue.

Point = zeros(Ktype+1,2);
if PCType == 1   % Square lattice
    Point(1,:)=[0,0];   %Gamma Point
    Point(2,:)=[b1(1)/2,0]; %X Point
    Point(3,:)=[(b1(1)+b2(1))/2,(b1(2)+b2(2))/2];   %M point
    Point(4,:)=[0,0];   %Gama Point
end
if PCType == 2  % Triangular lattice
    Point(1,:)=[0,0];   %Gamma Point
    Point(2,:)=[b1(1)/3,b1(2)]; %K Point
    Point(3,:)=[0,b1(2)];   %M point
    Point(4,:)=[0,0];   %Gama Point
end
K1 = [];
K2 = [];
K3 = [];
for ktnum=1:Ktype
    K1temp=linspace(Point(ktnum,1),Point(ktnum+1,1),Keach+1);
    K2temp=linspace(Point(ktnum,2),Point(ktnum+1,2),Keach+1);
    K1=[K1,K1temp(1:Keach)];
    K2=[K2,K2temp(1:Keach)];
end
K3 = kz_const*ones(size(K1));

eigval = [];    % Initialize the eigenvalue matrix.
for knum = 1:NumberK
    kx = K1(knum);
    ky = K2(knum);
    kz = 0;
    modulus_kG_m = sqrt((kx+Gx_m).^2 + (ky+Gy_m).^2 + kz^2);
    modulus_kG_n = sqrt((kx+Gx_n).^2 + (ky+Gy_n).^2 + kz^2);
    
    e1x_m = (ky+Gy_m)./modulus_kG_m;
    e1y_m = -(kx+Gx_m)./modulus_kG_m;
    % When kx+Gx_m, ky+Gy_m is zero, choose arbitrary e1_m
    e1x_m(isnan(e1x_m)) = 1/sqrt(2);
    e1y_m(isnan(e1y_m)) = 1/sqrt(2);
    e1z_m = zeros(size(Gx_m));
    
    e1x_n = (ky+Gy_n)./modulus_kG_n;
    e1y_n = -(kx+Gx_n)./modulus_kG_n;
    % When kx+Gx_n, ky+Gy_n is zero, choose arbitrary e1_n
    e1x_n(isnan(e1x_n)) = 1/sqrt(2);
    e1y_n(isnan(e1y_n)) = 1/sqrt(2);
    e1z_n = zeros(size(Gx_n));
    
    kz = K3(knum);
    modulus_kG_m = sqrt((kx+Gx_m).^2 + (ky+Gy_m).^2 + kz^2);
    modulus_kG_n = sqrt((kx+Gx_n).^2 + (ky+Gy_n).^2 + kz^2);
    
    e2x_m = (e1y_m*kz - e1z_m.*(ky+Gy_m))./modulus_kG_m;
    e2y_m = (e1z_m.*(kx+Gx_m) - e1x_m*kz)./modulus_kG_m;
    e2z_m = (e1x_m.*(ky+Gy_m) - e1y_m.*(kx+Gx_m))./modulus_kG_m;
    
    e2x_n = (e1y_n*kz - e1z_n.*(ky+Gy_n))./modulus_kG_n;
    e2y_n = (e1z_n.*(kx+Gx_n) - e1x_n*kz)./modulus_kG_n;
    e2z_n = (e1x_n.*(ky+Gy_n) - e1y_n.*(kx+Gx_n))./modulus_kG_n;
    
    KGmn_mat11 = modulus_kG_m.*modulus_kG_n.*(e2x_m.*e2x_n + e2y_m.*e2y_n + e2z_m.*e2z_n);
    H11 = KGmn_mat11.*ek_mat;
    KGmn_mat12 = -modulus_kG_m.*modulus_kG_n.*(e2x_m.*e1x_n + e2y_m.*e1y_n + e2z_m.*e1z_n);
    H12 = KGmn_mat12.*ek_mat;
    KGmn_mat21 = -modulus_kG_m.*modulus_kG_n.*(e1x_m.*e2x_n + e1y_m.*e2y_n + e1z_m.*e2z_n);
    H21 = KGmn_mat21.*ek_mat;
    KGmn_mat22 = modulus_kG_m.*modulus_kG_n.*(e1x_m.*e1x_n + e1y_m.*e1y_n + e1z_m.*e1z_n);
    H22 = KGmn_mat22.*ek_mat;
    
    H = [H11,H12;H21,H22];
    
    eigvalue = sort(eig(H));
    %eigval(:,knum) = eigvalue(1:NEIG);
    eigval=[eigval,eigvalue(1:NEIG)];   
end

eigval=[eigval,eigval(:,1)];
eigval = sqrt(eigval)*a*0.5/pi;
eigval=real(eigval); %Complete All the things

% Plot the figures
for m=1:Ktype
   D(m)=sqrt((Point(m+1,1)-Point(m,1))^2+(Point(m+1,2)-Point(m,2))^2);
   xtemp(m,:)=linspace(0,D(m),Keach+1);
end
x=xtemp(1,1:Keach);
Dtotal=0;
for m=2:Ktype
   Dtotal=Dtotal+D(m-1);
   x=[x,xtemp(m,1:Keach)+Dtotal];
end
x=[x,xtemp(Ktype,Keach+1)+Dtotal];
x=x/max(x);

MaxB=1.0;
x1=x(Keach+1);
x2=x(Keach*2+1);

figure;
clf;

h=plot(x,eigval,'b-',[x1 x1],[0 MaxB],'k:',[x2 x2],[0 MaxB],'k:');
set(h,'LineWidth',2.0);

axis([0 1 0 MaxB]);
%h=xlabel('in-plane wave vector k_/_/');
h=ylabel('Normalized frequency (\omegaa/2\pic)');
set(h,'FontSize',14);

if PCType == 1
    titletext=strcat('Square lattice (R=',num2str(R),',kz=',num2str(kz_const),',ea=',num2str(ea),',eb=',num2str(eb),')');
    text(x(1)-0.02,-0.03, '\Gamma','FontSize',14);
    text(x1-0.02,-0.03, 'X','FontSize',14)
    text(x2-0.02,-0.03, 'M','FontSize',14)
    text(x(Keach*Ktype+1)-0.02,-0.03, '\Gamma','FontSize',14)    
end
if PCType == 2
    titletext=strcat('Triangular lattice (R=',num2str(R),',kz=',num2str(kz_const),',ea=',num2str(ea),',eb=',num2str(eb),')');
    text(x(1)-0.02,-0.03, '\Gamma','FontSize',14)
    text(x1-0.02,-0.03, 'K','FontSize',14)
    text(x2-0.02,-0.03, 'M','FontSize',14)
    text(x(Keach*Ktype+1)-0.02,-0.03, '\Gamma','FontSize',14)    
end

h=title(titletext);
set(h,'FontSize',14);
set(gca,'xtick',[]);



%tmpx = 1:NumberK;
%figure,plot(tmpx,eigval*2*pi,'linewidth',2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% counter=1;
% for ii=1:size(k0,1),
%    k=k0(ii,:);
%    K(:,1)=k(1)+G(:,1);
%    K(:,2)=k(2)+G(:,2);
%    K(:,3)=0;
%    %%%find the unit cell perpendicular to it in xy plane
%    %%%be sure to deal with the case of modulus(k)=0
%    %%%NaN in this case
%    %%%to save computing time, use absk=modulus(K)
%    absk=modulus(K);
%    e1=[K(:,2)./absk,-K(:,1)./absk,zeros(length(K),1)];
%    e1(isnan(e1))=1/sqrt(2); %%%when Kz is not zero,kx,ky is zero choose arbitrary e1
%    K(:,3)=k(3)+G(:,3);
%    %%%find the other perpendicular unit cell
%    %%%be sure to deal with the case of modulus(k)=0
%    %%%NaN in this case
%    absk=modulus(K);
%    e2=cross(e1, K);
%    e2=[e2(:,1)./absk,e2(:,2)./absk,e2(:,3)./absk];
%    e2(isnan(e2))=0;
%    %%%form the equation matrix, it should be 2N by 2N
%    %%%in this case we will have no TE and TM decoupling
%    M1=([absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]*[absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]').*eps2;
%    M2=([absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]*[absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]').*eps2;
%    M3=([absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]*[absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]').*eps2;
%    M4=([absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]*[absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]').*eps2;
%    M=[M1,M3;M2,M4];
% 
% 	E=sort(abs(eig(M)));
% 	freq(:,counter)=sqrt(abs(E(1:20))).*a./2./pi;
% 	display(sprintf('calculation of k=[%f,%f,%f] is finished',k(1),k(2),k(3)));
% 	counter=counter+1;
% end
% tmpx=1:length(k0);
% %plot(n,freq,'o'),hold on
% plot(tmpx,freq*2*pi,'linewidth',2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%title('Full band structure of a 2D triangular silica-air hole with a kz')

%axis([0,length(k0)+1,0,1.4])