pwm法计算二维光子晶体能带的通用程序

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('Tria

4、ngular lattice');endif (PCType=3) disp('Honeycomb lattice');end%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.a

5、1=a*1,0;if PCType=1 a2=a*0,1;endif PCType=2 a2=a*0.5,sqrt(3)/2;end%a1,a2 are the basic 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.

6、 the fraction of% the total volume occupied by the rods.MaxDimForG=12; % The max Potive 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

7、O O% O O O O O % | % |%initial the G matrixdisp('-')disp('Dielectric constant FT- BEGIN')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

8、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; 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)*

9、R;if PCType=1|PCType=2 %eliminate the division on zero in the calculatation of ek na=find(GR_mat=0); GR_mat(na)=1; ek_mat=ekc*besselj(1,GR_mat)./GR_mat; ek_mat(na)=ek0;endif PCType=3 %eliminate the division on zero in the calculatation of ek na=find(GR_mat=0); GR_mat(na)=1; ek_mat=cos(Gx_m-Gx_n).*a/

10、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 lattice, or rectangular lattice Point(1,:)=0,0; %Gama 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,:

11、)=0,0; %Gama Pointendif PCType=2|PCType=3 %Triangular lattice and Honeycomb lattice. Point(1,:)=0,0; %Gama Point Point(2,:)=b2(2)*sqrt(3)/6,b2(2)/2; %K Point Point(3,:)=0,b2(2)/2; % M point Point(4,:)=0,0; %Gama Pointend%These three are for the K vectors, for the different case.K1=;K2=;for ktnum=1:K

12、type 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);enddisp('Dielectric constant FT- END')disp('-')disp('Eigen value calculations- BEGIN')eigval=; %Initial the ei

13、gvalue matrix.for knum=1:NumberK disp(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 part KGmn_mat=(kx-Gx_m).*(kx-Gx_n)+(ky-Gy_m).*(ky-Gy_n); H=KGmn_mat.*ek_mat; else KGmn_mat=(kx-Gx_m).

14、*(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 eigenvalues eigvalue=sort(eig(H); eigval=eigval,eigvalue(1:NEIG); endeigval=eigval,eigval(:,1); eigval=sqrt(eigval)*a*0.5/pi;eigval=real(eigval); %Complete All the thingsdisp('Eigen va

15、lue calculation0s- END')disp('-')%Plot the figures%x=linspace(0,10,NumberK+1);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);endx=xtemp(1,1:Keach);Dtotal=0;for m=2:Ktype Dtotal=Dtotal+D(m-1); x=x,xtemp(m,1:Keach)+Dtotal;en

16、dx=x,xtemp(Ktype,Keach+1)+Dtotal;x=x/max(x);MaxB=0.8;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);if (TEorTM=0) legend('TE modes',4);else legend('TM modes',4);end axis(0 1 0 M

17、axB);h=ylabel('Normalized frequency (a/lambda)');set(h,'FontSize',14);if (PCType=1) titletext=strcat('Square Lattice (ea=',num2str(ea),', eb=',num2str(eb),', R=',num2str(R),')'); text(x(1)-0.02,-0.03, 'Gamma','FontSize',14) text(x1-0.02

18、,-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) endif (PCType=2) titletext=strcat('Triangular Lattice (ea=',num2str(ea),', eb=',num2str(eb),', R=',num2str(R),')'); 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','Fon