《pwm法计算二维光子晶体能带的通用程序》由会员分享，可在线阅读，更多相关《pwm法计算二维光子晶体能带的通用程序（6页珍藏版）》请在金锄头文库上搜索。
1、function PBGBand(ea,eb,R,PCType,Keach,TEorTM)%function PBGBand(ea,eb,R,PCType,Keach,TEorTM)%-%| This is a program to calculate the Photonic Bands of two |%| dimension Photonic Crystal with circular inclusions. |%| It calculates both TE and TM modes (E and H polarization) |%-%Parameters:%ea: The diel
2、ectric constant of the circular inclusions.%eb: The dielectric constant of background.%R: The radius of dielectric columns%PCType =1: Square lattice% =2: Triangular lattice% =3: Honeycomb%Keach: The number of k vectors in each wave vector branch.%TEorTM: =0: TE modes% =1: TM modes%-tic;TEorTM=1;PCTy
3、pe=2;Keach=6;R=0.3;ea=1;eb=10.5;disp(-)if (TEorTM=0)disp(Plane wave expansion method for PC bands: TE modes);else;disp(Plane wave expansion method for PC bands: TM modes);end disp(-)if (PCType=1)disp(Square lattice);endif (PCType=2)disp(Triangular lattice);endif (PCType=3)disp(Honeycomb lattice);end
4、%Control parametersKtype=3; % The number of band parts, such as X-M, T-X, .NumberK=Ktype*Keach; %The total number of K vector;NEIG=20; %NEIG: The number cut of the Eigvalue.%Initial parametersa=1; %Lattice constance.a1=a*1,0;if PCType=1a2=a*0,1;endif PCType=2a2=a*0.5,sqrt(3)/2;end%a1,a2 are the basi
5、c vectors of lacctice cell.ac=abs(a1(1)*a2(2)-a1(2)*a2(1); %ac: Area of lattice cell.b1=2*pi/ac*a2(2),-a2(1);b2=2*pi/ac*-a1(2),a1(1);%b1, b2 are vectors in reciprocal space.f=pi*R*R/ac;%f: The filling fraction, i.e. the fraction of% the total volume occupied by the rods.MaxDimForG=12; % The max Poti
6、ve Number of the reciprocal lattice, GDimForG=2*MaxDimForG+1; NPW=DimForG*DimForG; %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 % | % |%initial the G matrixdisp(-)disp(Dielectric constant FT- BEGIN)gtemp=-Max
7、DimForG: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; %calculate the Matrix coefficience.ek0=f/ea+(1-f)/eb;
8、ekc=(1/ea-1/eb)*f*2; %Calculate the ek matrix, the coefficence of Fourier transform of ek.GR_mat=sqrt(Gx_m-Gx_n).*(Gx_m-Gx_n)+(Gy_m-Gy_n).*(Gy_m-Gy_n)*R;if PCType=1|PCType=2%eliminate the division on zero in the calculatation of ekna=find(GR_mat=0);GR_mat(na)=1; ek_mat=ekc*besselj(1,GR_mat)./GR_mat;
9、 ek_mat(na)=ek0;endif PCType=3%eliminate the division on zero in the calculatation of ekna=find(GR_mat=0);GR_mat(na)=1; ek_mat=cos(Gx_m-Gx_n).*a/2+(Gy_m-Gy_n).*a*sqrt(3)/6).*ekc.*besselj(1,GR_mat)./GR_mat; ek_mat(na)=ek0;end%toc%tic%Calculated points:Point=zeros(Ktype+1,2);if PCType=1 %Square lattic
10、e, or rectangular lattice Point(1,:)=0,0; %Gama PointPoint(2,:)=b1(1)/2,0; %X PointPoint(3,:)=(b1(1)+b2(1)/2,(b1(2)+b2(2)/2; %M pointPoint(4,:)=0,0; %Gama Pointendif PCType=2|PCType=3 %Triangular lattice and Honeycomb lattice.Point(1,:)=0,0; %Gama PointPoint(2,:)=b2(2)*sqrt(3)/6,b2(2)/2; %K PointPoi
11、nt(3,:)=0,b2(2)/2; % M pointPoint(4,:)=0,0; %Gama Pointend%These three are for the K vectors, for the different case.K1=;K2=;for ktnum=1:KtypeK1temp=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);
12、enddisp(Dielectric constant FT- END)disp(-)disp(Eigen value calculations- BEGIN)eigval=; %Initial the eigvalue matrix.for knum=1:NumberKdisp(strcat(-K vector No.,num2str(knum),-,num2str(NumberK) kx=K1(knum);ky=K2(knum);%tic%Now begin to calculate the matrix H: if (TEorTM=0)%TE partKGmn_mat=(kx-Gx_m)
13、.*(kx-Gx_n)+(ky-Gy_m).*(ky-Gy_n);H=KGmn_mat.*ek_mat;elseKGmn_mat=(kx-Gx_m).*(kx-Gx_n)+(ky-Gy_m).*(ky-Gy_n);H=KGmn_mat.*ek_mat;%TM part %Place your codes below for Task 1.end %Find the eigenvalueseigvalue=sort(eig(H);eigval=eigval,eigvalue(1:NEIG); endeigval=eigval,eigval(:,1); eigval=sqrt(eigval)*a*
14、0.5/pi;eigval=real(eigval); %Complete All the thingsdisp(Eigen value calculation0s- END)disp(-)%Plot the figures%x=linspace(0,10,NumberK+1);for m=1:KtypeD(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);endx=xtemp(1,1:Keach);Dtotal=0;for m=2:KtypeDtotal=Dtotal+D(m-1);x=x,xtemp(m,1:Keach)+Dtotal;endx=x,xtemp(Ktype,
