一個二維的FDTD程序% 本程序實現(xiàn)2維TM波FDTD仿真% 此程序用PML設置吸收邊界條件% FDTD_2D_kongqi_PML% 僅含有Ez,Hx,Hy分量clear;clc;% 1.初始化T=200; % 迭代次數(shù) IE=100; % JE=100;npml=8; % PML的網(wǎng)格數(shù)量c0=3*108; % 波速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*109); % 真空介電常數(shù)epsilon=1; % 相對介電常數(shù)sigma=0; % 電導率spread=6; % 脈沖寬度t0=20; % 脈沖高度ic=IE/2; % 源的X位置jc=JE/2; % 源的Y位置for i=1E+1;for j=1:JE+1;dz(i,j)=0; % z方向電荷密度ez(i,j)=0; % z方向電場hx(i,j)=0; % x方向磁場hy(i,j)=0; % y方向磁場 ihx(i,j)=0;%ihy(i,j)=0;iz(i,j)=0; % z方向求和參量，頻域卷積轉化為時域求和end;end;for i=2E; % for j=2:JE;ga(i,j)=1;end;end; %PML參數(shù)的設置for i=1E; gi2(i)=1;gi3(i)=1;fi1(i)=0;fi2(i)=1.0;fi3(i)=1.0;endfor j=1:JE;gj2(j)=1;gj3(j)=1;fj1(j)=0;fj2(j)=1;fj3(j)=1;endfor i=1:npml+1; %設置PML層中的參數(shù)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);endfor i=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.迭代求解電場和磁場for t=1:T;for i=2E; % 為了使每個電場周圍都有磁場進行數(shù)組下標處理for j=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; % 電場循環(huán)結束pulse=sin(2*pi*f*t*dt); % 正弦波源dz(ic,jc)=dz(ic,jc)+pulse; % 軟源for i=1E; % 為了使每個電場周圍都有磁場進行數(shù)組下標處理for j=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; % 電荷密度循環(huán)結束for j=1:JE;ez(1,j)=0;ez(IE,j)=0;endfor i=1E;ez(i,1)=0;ez(i,JE)=0;end;for i=1E; % 為了使每個磁場周圍都有電場進行數(shù)組下標處理for j=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循環(huán)結束for i=1E-1; % 為了使每個磁場周圍都有電場進行數(shù)組下標處理for j=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循環(huán)結束end;end; 在Maxwell旋度方程的差分表示中，按照Yee氏的空間網(wǎng)格設置，將出現(xiàn)半空間步長，通過前一時刻的磁、電場值得到時刻的電、磁場值，并在每一時刻上，將此過程算遍整個空間中隨時間變化的電、磁場值的解，但在編程計算中不使用1/2的空間表示，而要通過一定的相互關系把它表達出來，在自制的C程序中采用數(shù)組來表示上面的表達式中各場值及系數(shù)，其表達式在程序中表示如下所示：ezk+1ij=CAij*ezkij+CBij*CD*(hyki+1j-hykij （6）+hxkij-hxkij+1)hxk+1ij+1=hxkij+1+CD*(ezkij-ezkij+1) （7）hyk+1i+1j=hyki+1j+CD*(ezki+1j-ezkij) （8） 由于各場值起始均賦零，時間步數(shù)從零開始，每一時間步均按上面的順序在整個模擬區(qū)計算一遍，這樣場的實際那關系和空間關系就完全被體現(xiàn)出來。首先我們看空間關系，比較（6）（8）和（3）（5）很容易看處，在（6）中的hyij實質上代表的是，即對于而言i實際上代表。比較（7）和（4）亦可看出，hxij+1實際上代表的是，即對于而言j1實際上代表。這樣用此處的hxij+1,hxij,hyi+1j和hyij一起代入（6）中計算下一時間步中的電場值ezij。就是說雖然程序中沒有出現(xiàn)網(wǎng)格的中間點，但這些點上的場值實際上被計算了出來。現(xiàn)在來考察時間關系，還是從磁場分量的計算開始，假設