% 本程序?qū)崿F(xiàn)2維TM波FDTD仿真



% 此程序用PML設(shè)置吸收邊界條件



% 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*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); % 時(shí)間間隔






epsz=1/(4*pi*9*10^9); % 真空介電常數(shù)



epsilon=1; % 相對介電常數(shù)



sigma=0; % 電導(dǎo)率






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方向求和參量，頻域卷積轉(zhuǎn)化為時(shí)域求和



end;



end;






for i=2E; % 



for j=2:JE;



ga(i,j)=1;



end;



end; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%PML參數(shù)的設(shè)置



for i=1E; 



gi2(i)=1;



gi3(i)=1;



fi1(i)=0;



fi2(i)=1.0;



fi3(i)=1.0;



end






for j=1:JE;



gj2(j)=1;



gj3(j)=1;



fj1(j)=0;



fj2(j)=1;



fj3(j)=1;



end






for i=1:npml+1; %設(shè)置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);



end






for 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; % 為了使每個(gè)電場周圍都有磁場進(jìn)行數(shù)組下標(biāo)處理



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)結(jié)束






pulse=sin(2*pi*f*t*dt); % 正弦波源



dz(ic,jc)=dz(ic,jc)+pulse; % 軟源






for i=1E; % 為了使每個(gè)電場周圍都有磁場進(jìn)行數(shù)組下標(biāo)處理



for j=1:JE;



ez(i,j)=ga(i,j)* dz(i,j); %反映煤質(zhì)的情況都是放到這里的



% iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ;



end;



end; % 電荷密度循環(huán)結(jié)束






for j=1:JE;



ez(1,j)=0;



ez(IE,j)=0;



end






for i=1E;



ez(i,1)=0;



ez(i,JE)=0;



end;






for i=1E; % 為了使每個(gè)磁場周圍都有電場進(jìn)行數(shù)組下標(biāo)處理



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)結(jié)束






for i=1E-1; % 為了使每個(gè)磁場周圍都有電場進(jìn)行數(shù)組下標(biāo)處理



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)結(jié)束



end;



end;