一个fdtd程序

一个fdtd程序以下是以一个方柱为例来仿真其散射的波形：程序如下：clear;clc;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%c=3*10^8;%波速f=3*10^9;%频率lamda=c/f;%波长k=2*pi/lamda;%波数epsz=1/(4*pi*9*10^9);%真空介电常数mu=4*pi*10.^(-7);%真空磁导率Z=120*pi;%真空波阻抗epsilon=1;%相对介电常数sigma=0;%电导率N=100;%网格数量L=800;%迭代次数ddx=lamda/20;%网格尺寸dt=ddx/(2*c);%时间间隔ia=N/4;%总场区域x左ib=3*N/4;%总场区域x右ja=ia;%总场区域x下jb=ib;%总场区域x上npml=N/8;%PML点数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%r=2*lamda;M=r/(2*ddx);fori=1:Nforj=1:Nga(i,j)=1/(epsilon+sigma*dt/epsz);%求和参量gb(i,j)=sigma*dt/epsz;%求和参量endendspread=8;%脉冲宽度t0=25;%脉冲高度is=N/2;%源的X位置js=N/2;%源的Y位置ez_inc=zeros(1,N);hx_inc=zeros(1,N);ez_low1=0;ez_low2=0;ez_low3=0;ez_low4=0;dz=zeros(N,N);%z方向电荷密度ez=zeros(N,N);%z方向电场iz=zeros(N,N);%z方向电场求和参量hx=zeros(N,N);%x方向磁场hy=zeros(N,N);%y方向磁场ihx=zeros(N,N);%x方向磁场参量ihy=zeros(N,N);%y方向磁场参量%%%%%%%%%%%%%%%%%%%%%%%%PML%%%%%%%%%%%%%%%%%%%fori=1:Ngi2(i)=1;gi3(i)=1;fi1(i)=0;fi2(i)=1;fi3(i)=1;endforj=1:Ngj2(j)=1;gj3(j)=1;fj1(j)=0;fj2(j)=1;fj3(j)=1;endfori=1:npml+1xnum=npml-i+1;xxn=xnum/npml;xn=0.33*(xxn^3);gi2(i)=1/(1+xn);gi2(N-i+1)=1/(1+xn);gi3(i)=(1-xn)/(1+xn);gi3(N-i+1)=(1-xn)/(1+xn);xxn=(xnum-0.5)/npml;xn=0.25*(xxn^3);fi1(i)=xn;fi1(N-i)=xn;fi2(i)=1/(1+xn);fi2(N-i)=1/(1+xn);fi3(i)=(1-xn)/(1+xn);fi3(N-i)=(1-xn)/(1+xn);endforj=1:npml+1xnum=npml-j+1;xxn=xnum/npml;xn=0.33*(xxn^3);gj2(j)=1/(1+xn);gj2(N-j+1)=1/(1+xn);gj3(j)=(1-xn)/(1+xn);gj3(N-j+1)=(1-xn)/(1+xn);xxn=(xnum-0.5)/npml;xn=0.25*(xxn^3);fj1(j)=xn;fj1(N-j)=xn;fj2(j)=1/(1+xn);fj2(N-j)=1/(1+xn);fj3(j)=(1-xn)/(1+xn);fj3(N-j)=(1-xn)/(1+xn);endforT=1:L%%%%%%%%%%%%%%TM波Y方向传播%%%%%%%%%%%%%%%%%%%%%forj=2:N-1ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));endez_inc(1)=ez_low2;ez_low2=ez_low1;ez_low1=ez_inc(2);ez_inc(N)=ez_low3;ez_low3=ez_low4;ez_low4=ez_inc(N-1);%%%%%%%%%%%%%%%%电荷密度dz%%%%%%%%%%%%%%%%%%%%%%%%%fori=2:N-1forj=2:Ndz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));endend%%%%%%%%%%%%%%%%%%脉冲的加入%%%%%%%%%%%%%%%%%%%source=exp((-0.5)*((t0-T)/spread).^2);ez_inc(5)=source;fori=ia:ibdz(i,ja)=dz(i,ja)+0.5*hx_inc(ja-1);dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb);end%%%%%%%%%%%%%%电场ez%%%%%%%%%%%%%%%%%%%%%%%%fori=1:Nforj=1:Nez(i,j)=ga(i,j)*(dz(i,j)-iz(i,j));iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j);endend%%%%%%%%%%%%%%%%%%%%边界电场%%%%%%%%%%%%%%%%%%%%%%%forj=1:Nez(1,j)=0;ez(N,j)=0;endfori=1:Nez(i,1)=0;ez(i,N)=0;end%%%%%%%%%%%%%%%%%%%%%%边界条件%%%%%%%%%%%%%%%%%%%%%%%%fori=N/2-M:N/2+M-1forj=N/2-M:N/2+M-1ez(i,j)=0;endend%%%%%%%%%%%%%%%%%%%%%%%%%磁场%%%%%%%%%%%%%%%%%%%%%%%%forj=1:N-1hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1));end%%%%%%%%%%%%%%%%%%%%%%%%磁场Hx%%%%%%%%%%%%%%%%%%%%%%fori=1:Nforj=1:N-1curl_e=ez(i,j)-ez(i,j+1);ihx(i,j)=ihx(i,j)+fi1(i)*curl_e;hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j));endend;fori=ia:ibhx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja);hx(i,jb)=hx(i,jb)-0.5*ez_inc(jb);end%%%%%%%%磁场Hy%%%%%%%%%%%%%%%%%%%fori=1:N-1forj=1:Ncurl_e=ez(i+1,j)-ez(i,j);ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j));endendforj=ja:jb;hy(ia-1,j)=hy(ia-1,j)-0.5*ez_inc(j);hy(ib,j)=hy(ib,j)+0.5*ez_inc(j);endmesh(ez)drawnow;end