该程序用来分析自由空间中的电磁波传输情况

%该程序用来分析自由空间中的电磁波传输情况

%该电磁波为TM 波

%激励源可以是高斯脉冲或者正弦波

%该仿真中加入了PML 边界

clear all

%设置网格多少

IE = 100; % x 方向网格多少

JE = 100; % y 方向网格多少

% %设置脉冲函数中心

ic = IE / 2;

jc = JE / 2;

%初始化设置

ddx = 0.01; %空间网格大小，x 方向和y 方向网格大小相同

dt = ddx / 6e8; %时域网格大小

epsz = 8.8e-12; %介电常数

pi = 3.13159; %常数PI

er=4; %介质柱相对介电常数

%初始化系统变量

dz = zeros(IE,JE); %电场通量D

hx = dz; %x 方向磁场强度

hy = dz; %y 方向磁场强度

ihy = dz; %用于中间变量ihy，是用来计算磁场强度的，是curl_e 的积分

ihx = dz; %用于中间变量ihy，是用来计算磁场强度的，是curl_e 的积分

ga = ones(IE,JE); %电场强度与D 之间的关系矩阵.

ez_inc=zeros(IE);

hx_inc=zeros(IE);

%确定介质柱

radius=0.2;

ia=30;%确定总场与散射场区域

ib=70;

ja=30;

jb=70;

for j=ja:jb

for i=ia:ib

xdist=ic-i;

ydist=jc-j;

dist=sqrt(xdist^2+ydist^2);

if dist<=radius

ga(i,j)=1./er;

end

end

end

%计算PML 参数，若fj1 和fi1 为0，而其他参数为1，则为自由空间中的情况

for i = 1:IE

gi2(i) = 1;

gi3(i) = 1;

fi1(i) = 0;

fi2(i) = 1;

fi3(i) = 1;

end;

for j = 1:JE

gj2(j) = 1;

gj3(j) = 1;

fj1(j) = 0;

fj2(j) = 1;

fj3(j) = 1;

end;

%输入PML Cell 个数，即PML 有多少个单元网格，在此，x 方向和y 方向上的PML 网格相同

npml = input('Please input the number of PML Cell: ');

%x 方向用*i*表示，一方向用*j*表示

%x 方向上的PML 参数设置

for i = 1:npml+1

xnum = npml -i + 1; %从npml 到0

xd = npml;

xxn = xnum / xd; %辅助变量xxn，从1 到0

xn = 0.33 * xxn^3; %成立方衰减

gi2(i) = 1 / ( 1 + xn );

gi2(IE-i+1) = 1 / ( 1 + xn );

gi3(i) = (1 - xn) / (1 + xn);

gi3(IE-i+1) = ( 1 - xn ) / ( 1 + xn );

xxn = ( xnum - 0.5 ) / xd;

xn = 0.25 * xxn^3;

fil(i) = xn;

fil(IE+1-i) = xn;

fi2(i) = 1 / ( 1 + xn );

fi2(IE+1-i) = 1 / ( 1 + xn );

fi3(i) = ( 1-xn ) / ( 1 + xn );

fi3(IE+1-i) = ( 1 - xn ) / ( 1 + xn );

end;

for j = 1:npml+1

xnum = npml - j + 1;

xd = npml;

xxn = xnum / xd;

xn = 0.33 * xxn^3;

gj2(j) = 1 / ( 1 + xn );

gj2(JE+1-j) = 1 / ( 1 + xn );

gj3(j) = ( 1 - xn ) / ( 1 + xn );

gj3(JE-j+1) = ( 1 - xn ) / ( 1 + xn );

xxn = ( xnum - 0.5 ) / xd;

xn = 0.25 * xxn^3;

fj1(j) = xn;

fj1(JE+1-j) = xn;

fj2(j) = 1 / ( 1 + xn );

fj2(JE+1-j) = 1 / ( 1 + xn );

fj3(j) = ( 1 - xn ) / ( 1 + xn );

fj3(JE+1-j) = ( 1 - xn ) / ( 1 + xn );

end;

%高斯脉冲变量设置

t0 = 10;

spread = 10;

T = 0;

%输入nsteps 必须为正整数

nsteps = input('Please input the number of nsteps ');

%设置循环次数，从常数可以得到空间和时间上的传输情况

%通过修改对应的部分数值即可得到书Fig3.2 图

%nsteps = 50

% nsteps = 20;

% nsteps = 30;

% nsteps = 40;

%主循环

for n = 1:nsteps

T = T+1;

%计算入射电场

for i=2:JE

ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));

end

%计算D

for j = 2:JE

for i = 2:IE

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;

% 激励元 在此采用正弦波，频率为1.5G，可选为高斯脉冲

% hard source

pulse = exp( -0.5 * (( t0-T ) / spread ) ^ 2 ) ;

dz(IE/2,JE/2) = pulse + dz(IE/2,JE/2);

ez_inc(50)=pulse;

%计算E

for j = 1:JE

for i = 1:IE

ez(i,j) = ga(i,j) * dz(i,j);

end;

end;

%Set the Ez edges to 0, as part of the PML

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;

%计算Hx

for j = 1:JE-1

for i = 1:IE

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;

%计算Hy

for j = 1:JE

for i = 1:IE -1

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;

%画图

figure(1);

surfc(1:IE,1:JE,ez)

axis([0 IE 0 JE -0.2 1.])

fmat=getframe(20)

end;