%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);
评论1