4242阅读
4回复

[转载]详细讨论一下这个二维的fdtd程序

UID ：3862

注册：2007-07-15
登录：2008-06-06
发帖：8
等级：旁观者

0楼发表于: 2007-10-24 19:20:06

clear;clc;

% 1.初始化

T=200; % 迭代次数
Nx=100;
Ny=100;
npml=8; %PML的网格数量

c0=3*10^8;
f=1.5*10^(9);
lambda=c0/f;

wl=10;
dx=lambda/wl; %取lamda/10的精度
dy=lambda/wl;
pi=3.14159;

dt=dx/(2*c0); %由于是二维情况同，取t=dx/2c0
t0=20*dt; %为高斯脉冲的仿真作准备

epsz=1/(4*pi*9*10^9);
epsilon=1;
sigma=0;

ic=Nx/4; % 源的X位置
jc=Ny/4; % 源的Y位置

for i=1:Nx+1;
for j=1:Ny+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=2:Nx;
for j=2:Ny;
ga(i,j)=1;
end;
end;
%PML参数的设置
for i=1:Nx;
gi2(i)=1; %这些数组代表是什么？
gi3(i)=1; %这些数组代表是什么？
fi1(i)=0; %这些数组代表是什么？
fi2(i)=1; %这些数组代表是什么？
fi3(i)=1; %这些数组代表是什么？
end

for j=1:Ny;
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; %这里的0.333并不是严格的计算，而是经验值
gi2(i)=1.0/(1+xn); %这些数组代表是什么？
gi2(Nx-1-i)=1/(1+xn); %这些数组代表是什么？
gi3(i)=(1-xn)/(1+xn); %这些数组代表是什么？
gi3(Nx-1-i)=(1-xn)/(1+xn); %这些数组代表是什么？

xn=0.25*((xnum-0.5)/npml)^3; %这里的0.25和0.333也是一个道理
fi1(i)=xn; %这些数组代表是什么？
fi1(Nx-2-i)=xn;
fi2(i)=1.0/(1+xn);
fi2(Nx-2-i)=1/(1+xn);
fi3(i)=(1-xn)/(1+xn);
fi3(Nx-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(Ny-1-i)=1/(1+xn);
gj3(i)=(1-xn)/(1+xn);
gj3(Ny-1-i)=(1-xn)/(1+xn);

xn=0.25*((xnum-0.5)/npml)^3;
fj1(i)=xn;
fj1(Ny-2-i)=xn;

fj2(i)=1.0/(1+xn);
fj2(Ny-2-i)=1/(1+xn);
fj3(i)=(1-xn)/(1+xn);
fj3(Ny-2-i)=(1-xn)/(1+xn);
end

%2.迭代求解电场和磁场

for t=1:T;
for i=2:Nx; % 为了使每个电场周围都有磁场进行数组下标处理
for j=2:Ny;
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); % 正弦波源
%pulse=exp(-4*pi*(t*dt-t0)^2/(2/f)^2);高斯脉冲情况
dz(ic,jc)=dz(ic,jc)+pulse; %软源，就是为了防止反射做的一个处理

for i=1:Nx; %为了使每个电场周围都有磁场进行数组下标处理
for j=1:Ny;
ez(i,j)=ga(i,j)* dz(i,j); %反映煤质的情况都是放到这里的
end;
end; % 电荷密度循环结束

for j=1:Ny;
ez(1,j)=0;
ez(Nx,j)=0;
end

for i=1:Nx;
ez(i,1)=0;
ez(i,Ny)=0;
end

for i=1:Nx; % 为使每个磁场周围都有电场进行数组下标处理
for j=1:Ny-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;
& ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册