mmi_fdtd2d_tm.zip_FDTD MMI_PML边界_fdtd波导_波导_点光源

%Two dimensional FDTD algorithm with PML BC and point source. close all; clc;clear; cc=299792458; eps0=8.854187817620391e-12; mu0= 1.256637061435917e-6; lam0=1;ds=lam0/40; Sclf=1/2; dt=Sclf*ds/cc; Lx=2*lam0;Ly=7*lam0; Nx=round(Lx/ds); Ny=round(Ly/ds); Nt=round(Ly*lam0/cc/dt*3.5+100); Nt2=Nt-100; Ez=zeros(Nx,Ny); Hx=zeros(Nx,Ny-1); Hy=zeros(Nx-1,Ny); Ezmax=Ez;Ezmin=Ez; % mur=1;epsr=1;sigma=0;sigmam=0; % Set the distribution of the MMI epsr=3.5^2;dy1=50;dy2=0;dx1=20;dx2=20;dwg=5; x0wg=round(Nx/2); epsrdist=ones(size(Ez)); epsrdist(dx1:Nx-dx2,dy1:Ny-dy2)=epsr; epsrdist(x0wg+(-dwg:dwg),1:dy1)=epsr; dpml=10; IHx=1:Nx;JHx=1:Ny-1; IHy=1:Nx-1;JHy=1:Ny; IEz=(dpml+2):(Nx-dpml-1);JEz=(dpml+2):(Ny-dpml-1); IEzx1=2:dpml+1;JEzx1=2:Ny-1;IEzx2=(Nx-dpml):(Nx-1);JEzx2=2:Ny-1; IEzy1=2:Nx-1;JEzy1=2:dpml+1;IEzy2=2:Nx-1;JEzy2=(Ny-dpml):(Ny-1); % PML parameters ------------------------------------------------- %The Only 3 parameters for PML npml=3.5;R0pml=1e-10;dpml=10; sig0E=-cc*eps0*log(R0pml)/2/dpml^(npml+1)/ds; sig0M=-cc*mu0 *log(R0pml)/2/dpml^(npml+1)/ds; %--------------------------------------------------------------- % Split fields to the x-(_x1), x+(_x2), y-(_y1) and y+(_y2) boundaries. Ezx_x1=zeros(dpml,Ny-2);Ezx_x2=Ezx_x1;Ezy_x1=Ezx_x1;Ezy_x2=Ezy_x1; Ezx_y1=zeros(Nx-2,dpml);Ezx_y2=Ezx_y1;Ezy_y1=Ezx_y1;Ezy_y2=Ezx_y1; % Sigma and sigmam to the x-, x+, y- and y+ boundaries sigmax_x1=zeros(dpml,Ny-2);sigmax_x2=zeros(dpml,Ny-2); sigmaxm_x1=zeros(dpml,Ny);sigmam_x2=sigmaxm_x1; sigmay_y1=zeros(Nx-2,dpml);sigmay_y2=sigmay_y1; sigmaym_y1=zeros(Nx,dpml);sigmaym_y2=sigmay_y1; % Setting the sigmas according to the (r/d)^n npml=npml+1; % Warning: npml=npml+1 lpml1=(dpml:-1:1)'-1;lpml2=(1:dpml)'-1; sigmax_x1=repmat(sig0E*((lpml1+1/2).^npml-(lpml1-1/2).^npml),[1,Ny-2]); sigmax_x2=repmat(sig0E*((lpml2+1/2).^npml-(lpml2-1/2).^npml),[1,Ny-2]); sigmax_x1(dpml,:)=repmat(sig0E*(0+1/2).^npml,[1,Ny-2]); sigmax_x2(1,:)=repmat(sig0E*(0+1/2).^npml,[1,Ny-2]); sigmay_x1=0*sigmax_x1;sigmay_x2=0*sigmax_x2; lpml1=(dpml:-1:1)'-1/2;lpml2=(1:dpml)'-1/2; sigmaxm_x1=repmat(sig0M*((lpml1+1/2).^npml-(lpml1-1/2).^npml),[1,Ny]); sigmaxm_x2=repmat(sig0M*((lpml2+1/2).^npml-(lpml2-1/2).^npml),[1,Ny]); sigmaym_x1=0*sigmaxm_x1;sigmaym_x2=0*sigmaxm_x2; lpml1=(dpml:-1:1)-1;lpml2=(1:dpml)-1; sigmay_y1=repmat(sig0E*((lpml1+1/2).^npml-(lpml1-1/2).^npml),[Nx-2,1]); sigmay_y2=repmat(sig0E*((lpml2+1/2).^npml-(lpml2-1/2).^npml),[Nx-2,1]); sigmay_y1(:,dpml)=repmat(sig0E*(0+1/2).^npml,[Nx-2,1]); sigmay_y2(:,1)=repmat(sig0E*(0+1/2).^npml,[Nx-2,1]); sigmax_y1=0*sigmay_y1;sigmax_y2=0*sigmay_y2; lpml1=(dpml:-1:1)-1/2;lpml2=(1:dpml)-1/2; sigmaym_y1=repmat(sig0M*((lpml1+1/2).^npml-(lpml1-1/2).^npml),[Nx,1]); sigmaym_y2=repmat(sig0M*((lpml2+1/2).^npml-(lpml2-1/2).^npml),[Nx,1]); sigmaxm_y1=0*sigmaym_y1;sigmaxm_y2=0*sigmaym_y2; sigmaHx=zeros(size(Hx)); sigmaHx(:,1:dpml)=sigmaym_y1; sigmaHx(:,Ny-dpml:Ny-1)=sigmaym_y2; CHx1=(2*mu0-sigmaHx*dt)./(2*mu0+sigmaHx*dt); CHx2=2*dt/ds./(2*mu0+sigmaHx*dt); sigmaHy=zeros(size(Hy)); sigmaHy(1:dpml,:)=sigmaxm_x1; sigmaHy(Nx-dpml:Nx-1,:)=sigmaxm_x2; CHy1=(2*mu0-sigmaHy*dt)./(2*mu0+sigmaHy*dt); CHy2=2*dt/ds./(2*mu0+sigmaHy*dt); CEz1=1; CEz2=dt/ds/eps0./epsrdist(IEz,JEz); CEz3=dt/ds/eps0./epsrdist(IEz,JEz); xx=2*eps0*epsrdist(IEzx1,JEzx1)+sigmax_x1*dt; CEzx_x1_1=(2*eps0*epsrdist(IEzx1,JEzx1)-sigmax_x1*dt)./xx;CEzx_x1_2=2*dt./xx/ds; xx=2*eps0*epsrdist(IEzx2,JEzx2)+sigmax_x2*dt; CEzx_x2_1=(2*eps0*epsrdist(IEzx2,JEzx2)-sigmax_x2*dt)./xx;CEzx_x2_2=2*dt./xx/ds; xx=2*eps0*epsrdist(IEzx1,JEzx1)+sigmay_x1*dt; CEzy_x1_1=(2*eps0*epsrdist(IEzx1,JEzx1)-sigmay_x1*dt)./xx;CEzy_x1_2=2*dt./xx/ds; xx=2*eps0*epsrdist(IEzx2,JEzx2)+sigmay_x2*dt; CEzy_x2_1=(2*eps0*epsrdist(IEzx2,JEzx2)-sigmay_x2*dt)./xx;CEzy_x2_2=2*dt./xx/ds; xx=2*eps0*epsrdist(IEzy1,JEzy1)+sigmax_y1*dt; CEzx_y1_1=(2*eps0*epsrdist(IEzy1,JEzy1)-sigmax_y1*dt)./xx;CEzx_y1_2=2*dt./xx/ds; xx=2*eps0*epsrdist(IEzy2,JEzy2)+sigmax_y2*dt; CEzx_y2_1=(2*eps0*epsrdist(IEzy2,JEzy2)-sigmax_y2*dt)./xx;CEzx_y2_2=2*dt./xx/ds; xx=2*eps0*epsrdist(IEzy1,JEzy1)+sigmay_y1*dt; CEzy_y1_1=(2*eps0*epsrdist(IEzy1,JEzy1)-sigmay_y1*dt)./xx;CEzy_y1_2=2*dt./xx/ds; xx=2*eps0*epsrdist(IEzy2,JEzy2)+sigmay_y2*dt; CEzy_y2_1=(2*eps0*epsrdist(IEzy2,JEzy2)-sigmay_y2*dt)./xx;CEzy_y2_2=2*dt./xx/ds; t=1:Nt;tao=round(Nt/3);omg=2*pi*1/tao; Src=sin(t*dt*2*pi*cc/lam0); %Harmonic source %Src=10-15*cos(1*omg*t)+6*cos(2*omg*t)-cos(3*omg*t);Src(tao:end)=0; %Src=15*omg*sin(1*omg*t)-6*omg*2*sin(2*omg*t)+3*omg*sin(3*omg*t);Src(tao:end)=0; %Isrc=round(Nx/2);Jsrc=round(Ny/2); Isrc=x0wg+(-dwg:dwg);Jsrc=15; Ezt=zeros(1,Nt); for t=1:Nt Hx(IHx,JHx)=CHx1.*Hx(IHx,JHx)-CHx2.*(Ez(IHx,JHx+1)-Ez(IHx,JHx)); Hy(IHy,JHy)=CHy1.*Hy(IHy,JHy)+CHy2.*(Ez(IHy+1,JHy)-Ez(IHy,JHy)); Ez(IEz,JEz)=CEz1.*Ez(IEz,JEz)+CEz2.*(Hy(IEz,JEz)-Hy(IEz-1,JEz))-... CEz3.*(Hx(IEz,JEz)-Hx(IEz,JEz-1)); Ezx_x1=CEzx_x1_1.*Ezx_x1+CEzx_x1_2.*(Hy(IEzx1,JEzx1)-Hy(IEzx1-1,JEzx1)); Ezy_x1=CEzy_x1_1.*Ezy_x1-CEzy_x1_2.*(Hx(IEzx1,JEzx1)-Hx(IEzx1,JEzx1-1)); Ezx_x2=CEzx_x2_1.*Ezx_x2+CEzx_x2_2.*(Hy(IEzx2,JEzx2)-Hy(IEzx2-1,JEzx2)); Ezy_x2=CEzy_x2_1.*Ezy_x2-CEzy_x2_2.*(Hx(IEzx2,JEzx2)-Hx(IEzx2,JEzx2-1)); Ezx_y1=CEzx_y1_1.*Ezx_y1+CEzx_y1_2.*(Hy(IEzy1,JEzy1)-Hy(IEzy1-1,JEzy1)); Ezy_y1=CEzy_y1_1.*Ezy_y1-CEzy_y1_2.*(Hx(IEzy1,JEzy1)-Hx(IEzy1,JEzy1-1)); Ezx_y2=CEzx_y2_1.*Ezx_y2+CEzx_y2_2.*(Hy(IEzy2,JEzy2)-Hy(IEzy2-1,JEzy2)); Ezy_y2=CEzy_y2_1.*Ezy_y2-CEzy_y2_2.*(Hx(IEzy2,JEzy2)-Hx(IEzy2,JEzy2-1)); Ez(IEzx1,JEzx1)=Ezx_x1+Ezy_x1; Ez(IEzx2,JEzx2)=Ezx_x2+Ezy_x2; Ez(IEzy1,JEzy1)=Ezx_y1+Ezy_y1; Ez(IEzy2,JEzy2)=Ezx_y2+Ezy_y2; Ez(IEzx1,JEzy1)=Ezx_x1(IEzx1-1,JEzy1-1)+Ezy_y1(IEzx1-1,JEzy1-1); Ez(IEzx2,JEzy1)=Ezx_x2(IEzx1-1,JEzy1-1)+Ezy_y1(IEzx2-1,JEzy1-1); Ez(IEzx1,JEzy2)=Ezx_x1(IEzx1-1,JEzy2-1)+Ezy_y2(IEzx1-1,JEzy1-1); Ez(IEzx2,JEzy2)=Ezx_x2(IEzx1-1,JEzy2-1)+Ezy_y2(IEzx2-1,JEzy1-1); %Ez(Isrc,Jsrc)=Ez(Isrc,Jsrc)+Src(t); Ez(Isrc,Jsrc)=Src(t); %Ezt(t)=Ez(Isrc+5,Jsrc+5); if t>Nt2 Ezmax=max(Ezmax,Ez); Ezmin=min(Ezmin,Ez); end if mod(t,23)==0 %mesh(Ez);zlim([-1 1]*0.5);title(t);pause(0.1); imagesc(Ez);axis image xy; title(t); hold on;contour(0.1*epsrdist,1,'k');hold off;pause(0.1); end end %Ezamp=(Ezmax-Ezmin)/2; %figure;imagesc(Ezamp);axis image xy; %hold on;contour(0.1*epsrdist,1,'k');hold off;pause(0.1);