five.rar_PML_fdtd_matlab_一维电磁场_电磁场

s=0.5; nd=20; na=10; m=4; nx=100; ny=100; nz=100; ex=zeros(nx,ny+1,nz+1); ey=zeros(nx+1,ny,nz+1); ez=zeros(nx+1,ny+1,nz); hx=zeros(nx-1,ny,nz); hy=zeros(nx,ny-1,nz); hz=zeros(nx,ny,nz-1); dx=zeros(nx,ny-1,nz-1); dy=zeros(nx-1,ny,nz-1); dz=zeros(nx-1,ny-1,nz); dx1=dx; dy1=dy; dz1=dz; bx=hx; by=hy; bz=hz; bx1=bx; by1=by; bz1=bz; xe=4/5*s*(m+1)*((1:na-1)/na).^m; xm=4/5*s*(m+1)*((.5:na-.5)/na).^m; xe=[fliplr(xe),xe]; xm=[fliplr(xm),xm]; ae=(2-xe)./(2+xe); be=2./(2+xe); am=(2-xm)./(2+xm); bm=2./(2+xm); hfig=imagesc(0,0,ez(:,:,nx/2+1),0.001*[-1,1]); axis equal tight; colorbar; for n=1:1000 bx(:,:,na+1:end-na)=bx(:,:,na+1:end-na)-s*(diff(ez(2:end-1,:,na+1:end-na),1,2)-diff(ey(2:end-1,:,na+1:end-na),1,3)); bx(:,:,[1:na,end-na+1:end])=permute(am,[3,1,2]).*bx(:,:,[1:na,end-na+1:end])-s*permute(bm,[3,1,2]).*cat(3,diff(ez(2:end-1,:,1:na),1,2)-diff(ey(2:end-1,:,1:na+1),1,3),diff(ez(2:end-1,:,end-na+1:end),1,2)-diff(ey(2:end-1,:,end-na:end),1,3)); bx1(na:end-na+1,:,:)=bx(na:end-na+1,:,:)-bx1(na:end-na+1,:,:); bx1([1:na-1,end-na+2:end],:,:)=permute(1+xe/2,[2,3,1]).*bx([1:na-1,end-na+2:end],:,:)-permute(1-xe/2,[2,3,1]).*bx1([1:na-1,end-na+2:end],:,:); hx(:,na+1:end-na,:)=hx(:,na+1:end-na,:)+bx1(:,na+1:end-na,:); hx(:,[1:na,end-na+1:end],:)=permute(am,[1,2,3]).*hx(:,[1:na,end-na+1:end],:)+permute(bm,[1,2,3]).*bx1(:,[1:na,end-na+1:end],:); bx1=bx; by(na+1:end-na,:,:)=by(na+1:end-na,:,:)-s*(diff(ex(na+1:end-na,2:end-1,:),1,3)-diff(ez(na+1:end-na,2:end-1,:),1,1)); by([1:na,end-na+1:end],:,:)=permute(am,[2,3,1]).*by([1:na,end-na+1:end],:,:)-s*permute(bm,[2,3,1]).*cat(1,diff(ex(1:na,2:end-1,:),1,3)-diff(ez(1:na+1,2:end-1,:),1,1),diff(ex(end-na+1:end,2:end-1,:),1,3)-diff(ez(end-na:end,2:end-1,:),1,1)); by1(:,na:end-na+1,:)=by(:,na:end-na+1,:)-by1(:,na:end-na+1,:); by1(:,[1:na-1,end-na+2:end],:)=permute(1+xe/2,[1,2,3]).*by(:,[1:na-1,end-na+2:end],:)-permute(1-xe/2,[1,2,3]).*by1(:,[1:na-1,end-na+2:end],:); hy(:,:,na+1:end-na)=hy(:,:,na+1:end-na)+by1(:,:,na+1:end-na); hy(:,:,[1:na,end-na+1:end])=permute(am,[3,1,2]).*hy(:,:,[1:na,end-na+1:end])+permute(bm,[3,1,2]).*by1(:,:,[1:na,end-na+1:end]); by1=by; bz(:,na+1:end-na,:)=bz(:,na+1:end-na,:)-s*(diff(ey(:,na+1:end-na,2:end-1),1,1)-diff(ex(:,na+1:end-na,2:end-1),1,2)); bz(:,[1:na,end-na+1:end],:)=permute(am,[1,2,3]).*bz(:,[1:na,end-na+1:end],:)-s*permute(bm,[1,2,3]).*cat(2,diff(ey(:,1:na,2:end-1),1,1)-diff(ex(:,1:na+1,2:end-1),1,2),diff(ey(:,end-na+1:end,2:end-1),1,1)-diff(ex(:,end-na:end,2:end-1),1,2)); bz1(:,:,na:end-na+1)=bz(:,:,na:end-na+1)-bz1(:,:,na:end-na+1); bz1(:,:,[1:na-1,end-na+2:end])=permute(1+xe/2,[3,1,2]).*bz(:,:,[1:na-1,end-na+2:end])-permute(1-xe/2,[3,1,2]).*bz1(:,:,[1:na-1,end-na+2:end]); hz(na+1:end-na,:,:)=hz(na+1:end-na,:,:)+bz1(na+1:end-na,:,:); hz([1:na,end-na+1:end],:,:)=permute(am,[2,3,1]).*hz([1:na,end-na+1:end],:,:)+permute(bm,[2,3,1]).*bz1([1:na,end-na+1:end],:,:); bz1=bz; dx(:,:,na:end-na+1)=dx(:,:,na:end-na+1)+s*(diff(hz(:,:,na:end-na+1),1,2)-diff(hy(:,:,na:end-na+1),1,3)); dx(:,:,[1:na-1,end-na+2:end])=permute(ae,[3,1,2]).*dx(:,:,[1:na-1,end-na+2:end])+s*permute(be,[3,1,2]).*cat(3,diff(hz(:,:,1:na-1),1,2)-diff(hy(:,:,1:na),1,3),diff(hz(:,:,end-na+2:end),1,2)-diff(hy(:,:,end-na+1:end),1,3)); dx1(na+1:end-na,:,:)=dx(na+1:end-na,:,:)-dx1(na+1:end-na,:,:); dx1([1:na,end-na+1:end],:,:)=permute(1+xm/2,[2,3,1]).*dx([1:na,end-na+1:end],:,:)-permute(1-xm/2,[2,3,1]).*dx1([1:na,end-na+1:end],:,:); ex(:,na+1:end-na,2:end-1)=ex(:,na+1:end-na,2:end-1)+dx1(:,na:end-na+1,:); ex(:,[2:na,end-na+1:end-1],2:end-1)=permute(ae,[1,2,3]).*ex(:,[2:na,end-na+1:end-1],2:end-1)+permute(be,[1,2,3]).*dx1(:,[1:na-1,end-na+2:end],:); dx1=dx; dy(na:end-na+1,:,:)=dy(na:end-na+1,:,:)+s*(diff(hx(na:end-na+1,:,:),1,3)-diff(hz(na:end-na+1,:,:),1,1)); dy([1:na-1,end-na+2:end],:,:)=permute(ae,[2,3,1]).*dy([1:na-1,end-na+2:end],:,:)+s*permute(be,[2,3,1]).*cat(1,diff(hx(1:na-1,:,:),1,3)-diff(hz(1:na,:,:),1,1),diff(hx(end-na+2:end,:,:),1,3)-diff(hz(end-na+1:end,:,:),1,1)); dy1(:,na+1:end-na,:)=dy(:,na+1:end-na,:)-dy1(:,na+1:end-na,:); dy1(:,[1:na,end-na+1:end],:)=permute(1+xm/2,[1,2,3]).*dy(:,[1:na,end-na+1:end],:)-permute(1-xm/2,[1,2,3]).*dy1(:,[1:na,end-na+1:end],:); ey(2:end-1,:,na+1:end-na)=ey(2:end-1,:,na+1:end-na)+dy1(:,:,na:end-na+1); ey(2:end-1,:,[2:na,end-na+1:end-1])=permute(ae,[3,1,2]).*ey(2:end-1,:,[2:na,end-na+1:end-1])+permute(be,[3,1,2]).*dy1(:,:,[1:na-1,end-na+2:end]); dy1=dy; dz(nx/2,ny/2,nz/2+1)=dz(nx/2,ny/2,nz/2+1)+sin(2*pi*n*s/nd); dz(:,na:end-na+1,:)=dz(:,na:end-na+1,:)+s*(diff(hy(:,na:end-na+1,:),1,1)-diff(hx(:,na:end-na+1,:),1,2)); dz(:,[1:na-1,end-na+2:end],:)=permute(ae,[1,2,3]).*dz(:,[1:na-1,end-na+2:end],:)+s*permute(be,[1,2,3]).*cat(2,diff(hy(:,1:na-1,:),1,1)-diff(hx(:,1:na,:),1,2),diff(hy(:,end-na+2:end,:),1,1)-diff(hx(:,end-na+1:end,:),1,2)); dz1(:,:,na+1:end-na)=dz(:,:,na+1:end-na)-dz1(:,:,na+1:end-na); dz1(:,:,[1:na,end-na+1:end])=permute(1+xm/2,[3,1,2]).*dz(:,:,[1:na,end-na+1:end])-permute(1-xm/2,[3,1,2]).*dz1(:,:,[1:na,end-na+1:end]); ez(na+1:end-na,2:end-1,:)=ez(na+1:end-na,2:end-1,:)+dz1(na:end-na+1,:,:); ez([2:na,end-na+1:end-1],2:end-1,:)=permute(ae,[2,3,1]).*ez([2:na,end-na+1:end-1],2:end-1,:)+permute(be,[2,3,1]).*dz1([1:na-1,end-na+2:end],:,:); dz1=dz; set(hfig,'cdata',ez(:,:,nz/2+1)); drawnow; end saveas(gcf,[mfilename,'.png']);