2200阅读
0回复

[讨论]这个程序中ga(i,j)到底和电导率和介电常数是什么关系？

UID ：9226

注册：2008-03-11
登录：2008-03-12
发帖：1
等级：旁观者

0楼发表于: 2008-03-11 16:54:37

clear;clc;

%%%%%%%%%%%%%%% 1.初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

T=300; % 迭代次数
Nx=116;
Ny=116;
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/2; % 源的X位置
jc=Ny/2; % 源的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;

ga(8:Nx-8,8:60)=0.1;

%%%%%%%%%%%%%%%%%%%%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,108)=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;
end; % 磁场HX循环结束

for i=1:Nx-1; % 为了使每个磁 ..

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

发帖回复