complgxity.rar_源码

%This program calculate the wave guide by a photonic crystal. %For two dimension case. %Only for TM case, in this version. %Trial Version clear; %tic %Initial parameters and other things. W=0.36; %Normalized frequency %The following parameters are control parameters. WaveGuide=1; %If this program is for Wave Guide? If so, please specify 1. IsMovie=0; %If you want to play movie, please use 1. IsFigure=1; %If it will plot the figures? If so, please specify 1. WantToSeeEp=1; %Do you want to see the distrubution of Ep? If so, please specify 1. %End of defining control parameters MLatx=11; %How many Lattice cell in x direction. MLaty=11; %How many Lattice cell in y direction. NMlat=21; %The gird number in each Lattice Cell. %SHOULD BE ODD INTEGER!!!! if mod(NMlat,2)==0 NMlat=NMlat+1; end %Force it to be a odd integer! NTx=MLatx*NMlat+1; %It is the number of the Grid along x axis. NTy=MLaty*NMlat+1; %It is the number of the Grid along y axis. if WaveGuide==1 Nrow=4; %The row number of coloumns between the PML boundary and the waveguide. end NPML=12; %How many PML layers will be used in our computation. NTimeSteps=2500; %Total number of Time Steps Meach=20; %Define the interval for plot figures if IsFigure==1. %This also works for saving intervals. R=0.2; %The radius of dielectric columns ea=11.4; %The dielectric constant of these columns. Zmax=0.6; %The maximum value for z axis when plotting figures. Colormax=0.6; %The maximum value for colormap when plotting figures. %Some constants mu0=4*pi*1.0e-7; %Epsilon Zero, if using Gauss Unit, it equals to 1. e0=8.85*1e-12; %Mu Zero, if using Gauss Unit, it equals to 1. c=1/sqrt(mu0*e0); %The light speed. factor=mu0/e0; %The factor between conductivity and permeability. %Permeability=Conductivity*factor, in PML. a=1;%e-6; %The lattice constant. W=W*(2*pi*c/a); %frequency Dx=a/NMlat; %Delta x. Dy=Dx; %Delta y. Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy))/c; %Time interval %tic %In the following partm we define the dielectric constants: Ep=ones(NTx-1,NTy-1)*e0; Ep_cell=ones(NMlat,NMlat)*e0; x=-(NMlat-1)/2*Dx:Dx:(NMlat-1)/2*Dx; [X,Y]=meshgrid(x); X=X'; Y=Y'; flag=find(sqrt(X.^2+Y.^2)<R); Ep_cell(flag)=e0*ea; Ep=repmat(Ep_cell,MLatx,MLaty); if WaveGuide==1 eb=1; %The dielectric constant in the waveguide Ep(1:(MLatx-Nrow)*NMlat,Nrow*NMlat+1:(Nrow+1)*NMlat)=e0*eb; Ep((MLatx-Nrow-1)*NMlat+1:(MLatx-Nrow)*NMlat,Nrow*NMlat+1:MLaty*NMlat)=e0*eb; Ep((MLatx-Nrow-1)*NMlat+1:(MLatx-Nrow)*NMlat,Nrow*NMlat+1:(Nrow+1)*NMlat)=Ep_cell; Ep((MLatx-Nrow-2)*NMlat+1:(MLatx-Nrow-1)*NMlat,(Nrow+1)*NMlat+1:(Nrow+2)*NMlat)=e0*eb; end %toc if WantToSeeEp==1 x=0:Dx:(NMlat*MLatx-1)*Dx; x=x-(NMlat*MLatx-1)*Dx/2; y=0:Dy:(NMlat*MLaty-1)*Dy; y=y-(NMlat*MLaty-1)*Dy/2; [X,Y]=meshgrid(x,y); X=X'; Y=Y'; surf(X,Y,Ep/e0); shading interp; view(0,90); axis([min(x), max(x),min(y), max(y)]) axis off; disp('Press any key to continue...'); pause end %End of defining the Ep. if IsFigure==1 %Define the X Y coordinate for figures x=0:Dx:(NMlat*MLatx-1+2*NPML)*Dx; x=x-(NMlat*MLatx-1+2*NPML)*Dx/2; y=0:Dy:(NMlat*MLaty-1+2*NPML)*Dy; y=y-(NMlat*MLaty-1+2*NPML)*Dy/2; [X,Y]=meshgrid(x,y); X=X'; Y=Y'; end %Define the Ez, Hx, Hy, which are in the inside region. Ez=zeros(NTx-1,NTy-1); Hx=zeros(NTx-1,NTy); Hy=zeros(NTx,NTy-1); %Parameters about PML: n=4; %The order of the polynomial that decribes the conductivity profile. R=1e-10; Delta=NPML*Dx; SigmaMax=-(n+1)*e0*c*log(R)/(Delta*2); NUM=NPML*2:-1:1; Sigmax=SigmaMax*((NUM*Dx/2+Dx/2).^(n+1)-(NUM*Dx/2-Dx/2).^(n+1))/(Delta^n*Dx*(n+1)); Sigmay=Sigmax; SigmaBound=SigmaMax*(Dx/2).^(n+1)/(Delta^n*Dx*(n+1)); %Sigmax=SigmaMax*(NUM/6).^(n+1); %Another way, the definition. EzxPML1=zeros(NPML,NPML); EzyPML1=zeros(NPML,NPML); HxPML1=zeros(NPML,NPML); HyPML1=zeros(NPML,NPML); %Zone 1 Sigmax_z1=repmat(Sigmax(2:2:NPML*2)',1,NPML); Sigmax_x1=repmat(Sigmax(2:2:NPML*2)',1,NPML); Sigmax_y1=repmat(Sigmax(1:2:NPML*2-1)',1,NPML); Sigmay_z1=fliplr(repmat(Sigmax(2:2:NPML*2),NPML,1)); Sigmay_x1=fliplr(repmat(Sigmax(1:2:NPML*2-1),NPML,1)); Sigmay_y1=fliplr(repmat(Sigmax(2:2:NPML*2),NPML,1)); %Zone 1 EzxPML2=zeros(NPML,NPML); EzyPML2=zeros(NPML,NPML); HxPML2=zeros(NPML,NPML); HyPML2=zeros(NPML,NPML); %Zone 2 Sigmax_z2=flipud(Sigmax_z1); Sigmax_x2=flipud(Sigmax_x1); Sigmax_y2=flipud(Sigmax_y1); Sigmay_z2=Sigmay_z1; Sigmay_x2=Sigmay_x1; Sigmay_y2=Sigmay_y1; %Zone 2 EzxPML3=zeros(NPML,NPML); EzyPML3=zeros(NPML,NPML); HxPML3=zeros(NPML,NPML); HyPML3=zeros(NPML,NPML); %Zone 3 Sigmax_z3=Sigmax_z1; Sigmax_x3=Sigmax_x1; Sigmax_y3=Sigmax_y1; Sigmay_z3=fliplr(Sigmay_z1); Sigmay_x3=fliplr(Sigmay_x1); Sigmay_y3=fliplr(Sigmay_y1); %Zone 3 EzxPML4=zeros(NPML,NPML); EzyPML4=zeros(NPML,NPML); HxPML4=zeros(NPML,NPML); HyPML4=zeros(NPML,NPML); %Zone 4 Sigmax_z4=flipud(Sigmax_z1); Sigmax_x4=flipud(Sigmax_x1); Sigmax_y4=flipud(Sigmax_y1); Sigmay_z4=fliplr(Sigmay_z1); Sigmay_x4=fliplr(Sigmay_x1); Sigmay_y4=fliplr(Sigmay_y1); %Zone 4 EzxPMLA=zeros(NTx-1,NPML); EzyPMLA=zeros(NTx-1,NPML); HxPMLA=zeros(NTx-1,NPML); HyPMLA=zeros(NTx,NPML); %Zone A Sigmay_zA=repmat(Sigmay_z1(1,:),NTx-1,1); Sigmay_xA=repmat(Sigmay_x1(1,:),NTx-1,1); Sigmay_yA=repmat(Sigmay_y1(1,:),NTx,1); %Zone A EzxPMLB=zeros(NTx-1,NPML); EzyPMLB=zeros(NTx-1,NPML); HxPMLB=zeros(NTx-1,NPML); HyPMLB=zeros(NTx,NPML); %Zone B Sigmay_zB=fliplr(Sigmay_zA); Sigmay_xB=fliplr(Sigmay_xA); Sigmay_yB=fliplr(Sigmay_yA); %Zone B EzxPMLC=zeros(NPML,NTy-1); EzyPMLC=zeros(NPML,NTy-1); HxPMLC=zeros(NPML,NTy); HyPMLC=zeros(NPML,NTy-1);%Zone C Sigmax_zC=repmat(flipud(Sigmay_z1(1,:)'),1,NTy-1); Sigmax_xC=repmat(flipud(Sigmay_x1(1,:)'),1,NTy); Sigmax_yC=repmat(flipud(Sigmay_y1(1,:)'),1,NTy-1); %Zone C EzxPMLD=zeros(NPML,NTy-1); EzyPMLD=zeros(NPML,NTy-1); HxPMLD=zeros(NPML,NTy); HyPMLD=zeros(NPML,NTy-1);%Zone D Sigmax_zD=flipud(Sigmax_zC); Sigmax_xD=flipud(Sigmax_xC); Sigmax_yD=flipud(Sigmax_yC); %Zone D if IsMovie==1 Movie=moviein(NTimeSteps/Meach+1); Mnum=1; EzAll=[EzxPML3+EzyPML3,EzxPMLC+EzyPMLC,EzxPML1+EzyPML1; EzxPMLB+EzyPMLB,Ez,EzxPMLA+EzyPMLA; EzxPML4+EzyPML4,EzxPMLD+EzyPMLD,EzxPML2+EzyPML2]; HxAll=[HxPML3,HxPMLC,HxPML1;HxPMLB,Hx,HxPMLA;HxPML4,HxPMLD,HxPML2]; HyAll=[HyPML3,HyPMLC,HyPML1;HyPMLB,Hy,HyPMLA;HyPML4,HyPMLD,HyPML2]; figure(1) surf(X,Y,real(EzAll)) shading interp; axis([-0.5*a*MLatx-Dx*NPML 0.5*a*MLatx+Dx*NPML -0.5*a*MLaty-Dx*NPML 0.5*a*MLaty+Dx*NPML -Zmax Zmax]) zlabel('Ez'); Movie(:,1)=getframe; end if WaveGuide==1 Npy=NMlat*Nrow+round((NMlat+1)/2); Npx=NMlat*Nrow+round((NMlat+1)/2); else Npx=round(NTx/2); Npy=round(NTy/2); end expboundary=exp(-SigmaBound*Dt/e0); tic for m=1:NTimeSteps %tic %The following part is for the inside region. %H components. Hx(:,2:NTy-1)=Hx(:,2:NTy-1)-Dt*(Ez(:,2:NTy-1)-Ez(:,1:NTy-2))/(Dy*mu0); Hy(2:NTx-1,:)=Hy(2:NTx-1,:)+Dt*(Ez(2:NTx-1,:)-Ez(1:NTx-2,:))/(Dx*mu0); Hx(:,NTy)=expboundary*Hx(:,NTy)-(1-expboundary)*... (EzxPMLA(:,1)+EzyPMLA(:,1)-Ez(:,NTy-1))/(SigmaBound*factor*Dy); %Boundary A Hx(:,1)=expboundary*Hx(:,1)-(1-expboundary)*... (Ez(:,1)-EzxPMLB(:,NPML)-EzyPMLB(:,NPML))/(SigmaBound*factor*Dy); %Boundary B Hy(1,:)=expboundary*Hy(1,:)+(1-expboundary)*... (Ez(1,:)-EzxPMLC(NPML,:)-EzyPMLC(NPML,:))/(SigmaBound*factor*Dx); %Boundary C Hy(NTx,:)=expboundary*Hy(NTx,:)+(1-expboundary)*... (EzxPMLD(1,:)+EzyPMLD(1,:)-Ez(NTx-1,:))/(SigmaBound*factor*Dx); %Boundary D %The following part is for ZONE A. H components. HxPMLA(:,1:NPML-1)=exp(-Sigmay_xA(:,1:NPML-1)*Dt/e0).*HxPMLA(:,1