程序2：用FDTD方法分析二維TM波在正弦波源激勵情況下，正方形PML吸收邊界的性能
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位置

for i=1:IE+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;
end;
end;

for i=2:IE; % 
for j=2:JE;
ga(i,j)=1;
end;
end; 

%PML參數的設置
for i=1:IE; 
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; %設置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);
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 n=1:T;
for i=2:IE; 
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; % 電場循環結束
pulse=sin(2*pi*f*n*dt); % 正弦波源
dz(ic,jc)=dz(ic,jc)+pulse; 
for i=1:IE; 
for j=1:JE;
ez(i,j)=ga(i,j)* dz(i,j); %反映媒質的情況都是放到這里的
end;
end; % 電荷密度循環結束

for j=1:JE;
ez(1,j)=0;
ez(IE,j)=0;
end

for i=1:IE;
ez(i,1)=0;
ez(i,JE)=0;
end

for i=1:IE; 
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循環結束

for i=1:IE-1;
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循環結束
end;
mesh(ez);
drawnow;