一个二维的FDTD程序

1一个二维的FDTD程序%本程序实现2维TM波FDTD仿真%此程序用PML设置吸收边界条件%FDTD_2D_kongqi_PML%仅含有Ez,Hx,Hy分量clear;clc;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1.初始化T=200;%迭代次数IE=100;%JE=100;npml=8;%PML的网格数量c0=3*10^8;%波速f=1.5*10^(9);%频率lambda=c0/f;%波长wl=10;dx=lambda/wl;dy=lambda/wl;pi=3.14159;dt=dx/(2*c0);%时间间隔epsz=1/(4*pi*9*10^9);%真空介电常数epsilon=1;%相对介电常数sigma=0;%电导率spread=6;%脉冲宽度t0=20;%脉冲高度ic=IE/2;%源的X位置jc=JE/2;%源的Y位置fori=1E+1;forj=1:JE+1;dz(i,j)=0;%z方向电荷密度ez(i,j)=0;%z方向电场hx(i,j)=0;%x方向磁场hy(i,j)=0;%y方向磁场2ihx(i,j)=0;%ihy(i,j)=0;iz(i,j)=0;%z方向求和参量，频域卷积转化为时域求和end;end;fori=2E;%forj=2:JE;ga(i,j)=1;end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%PML参数的设置fori=1E;gi2(i)=1;gi3(i)=1;fi1(i)=0;fi2(i)=1.0;fi3(i)=1.0;endforj=1:JE;gj2(j)=1;gj3(j)=1;fj1(j)=0;fj2(j)=1;fj3(j)=1;endfori=1:npml+1;%设置PML层中的参数xnum=npml+1-i;xn=0.33*(xnum/npml)^3;gi2(i)=1.0/(1+xn);gi2(IE-1-i)=1/(1+xn);gi3(i)=(1-xn)/(1+xn);gi3(IE-1-i)=(1-xn)/(1+xn);xn=0.25*((xnum-0.5)/npml)^3;fi1(i)=xn;fi1(IE-2-i)=xn;fi2(i)=1.0/(1+xn);fi2(IE-2-i)=1/(1+xn);fi3(i)=(1-xn)/(1+xn);fi3(IE-2-i)=(1-xn)/(1+xn);end3fori=1:npml+1;xnum=npml+1-i;xn=0.33*(xnum/npml)^3;gj2(i)=1.0/(1+xn);gj2(JE-1-i)=1/(1+xn);gj3(i)=(1-xn)/(1+xn);gj3(JE-1-i)=(1-xn)/(1+xn);xn=0.25*((xnum-0.5)/npml)^3;fj1(i)=xn;fj1(JE-2-i)=xn;fj2(i)=1.0/(1+xn);fj2(JE-2-i)=1/(1+xn);fj3(i)=(1-xn)/(1+xn);fj3(JE-2-i)=(1-xn)/(1+xn);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2.迭代求解电场和磁场fort=1:T;fori=2E;%为了使每个电场周围都有磁场进行数组下标处理forj=2:JE;dz(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));end;end;%电场循环结束pulse=sin(2*pi*f*t*dt);%正弦波源dz(ic,jc)=dz(ic,jc)+pulse;%软源fori=1E;%为了使每个电场周围都有磁场进行数组下标处理forj=1:JE;ez(i,j)=ga(i,j)*dz(i,j);%反映煤质的情况都是放到这里的%iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j);end;end;%电荷密度循环结束forj=1:JE;ez(1,j)=0;ez(IE,j)=0;endfori=1E;ez(i,1)=0;4ez(i,JE)=0;end;fori=1E;%为了使每个磁场周围都有电场进行数组下标处理forj=1:JE-1;curl_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));end;end;%磁场HX循环结束fori=1E-1;%为了使每个磁场周围都有电场进行数组下标处理forj=1:JE;curl_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));end;end;%磁场HY循环结束end;end;5在Maxwell旋度方程的差分表示中，按照Yee氏的空间网格设置，将出现半空间步长，通过前一时刻的磁、电场值得到时刻的电、磁场值，并在每一时刻上，将此过程算遍整个空间中随时间变化的电、磁场值的解，但在编程计算中不使用1/2的空间表示，而要通过一定的相互关系把它表达出来，在自制的C程序中采用数组来表示上面的表达式中各场值及系数，其表达式在程序中表示如下所示：ez[k+1][i][j]=CA[i][j]*ez[k][i][j]+CB[i][j]*CD*(hy[k][i+1][j]-hy[k][i][j]（6）+hx[k][i][j]-hx[k][i][j+1])hx[k+1][i][j+1]=hx[k][i][j+1]+CD*(ez[k][i][j]-ez[k][i][j+1])（7）hy[k+1][i+1][j]=hy[k][i+1][j]+CD*(ez[k][i+1][j]-ez[k][i][j])（8）由于各场值起始均赋零，时间步数从零开始，每一时间步均按上面的顺序在整个模拟区计算一遍，这样场的实际那关系和空间关系就完全被体现出来。首先我们看空间关系，比较（6）－（8）和（3）－（5）很容易看处，在（6）中的hy[i][j]实质上代表的是),21(jiHy，即对于yH而言i实际上代表21i。比较（7）和（4）亦可看出，hx[i][j+1]实际上代表的是)21,(jiHx，即对于xH而言j＋1实际上代表21j。这样用此处的hx[i][j+1],hx[i][j],hy[i+1][j]和hy[i][j]一起代入（6）中计算下一时间步中的电场值ez[i][j]。就是说虽然程序中没有出现网格的中间点，但这些点上的场值实际上被计算了出来。现在来考察时间关系，还是从磁场分量的计算开始，假设现在是第k步，则在（6）－（8）6中hy[k][i][j]实际上表示的为),(21jiHky，hx[k][i][j]实际上表示的为),(21jiHkx，式中对于yH，xH而言，实际上k表示21k,首先由（7）和（8）推算出hx[k+1][i][j]，hy[k+1][i+1][j]，然后把这两式代入下一时间步中，由式（6）便可计算出ez[k+1][i][j]的值。这就完成了各场量的时间递推关系。