1379阅读
1回复

[资料共享]Electromagnetic Simulation using FDTD Sullivan书中的二维FDTD

人生若只如初见...

UID ：91870

注册：2012-04-11
登录：2012-06-03
发帖：46
等级：仿真新人

0楼发表于: 2012-05-06 21:55:18

                                                                                          %使用matlab程序实现
%Electromagnetic Simulation using FDTD Sullivan书中的二维FDTD
%by monkeyshen
%monkeyshenTech@qq.com
%2008年10月
clear all

%设置网格大小
IE = 60;
JE = 60;

%设置脉冲函数中心
ic = IE / 2;
jc = JE / 2;

%初始化设置
ddx = 0.01;         %空间网格大小
dt = ddx / 6e8;     %时域网格大小
epsz = 8.8e-12;     %介电常数
pi = 3.13159;

%初始化系统变量
dz = zeros(60,60);  %Electric Flux
ez = dz;            %Electric Field
hx = dz;            %Magnetic Field x vector
hy = dz;            %Magnetic Field y vector
ga = ones(60,60);

%for   j = 1:JE
%    for i = 1:IE
%        dz(i,j) = 0;
%        ez(i,j) = 0;
%        hx(i,j) = 0;
%        hy(i,j) = 0;
%        ga(i,j) = 1.0;
%    end;
%end;

%高斯脉冲变量设置
t0 = 20;
spread = 6;
T = 0;

%设置循环次数，从常数可以得到空间和时间上的传输情况
%通过修改对应的部分数值即可得到书Fig3.2图
nsteps = 50
% nsteps = 20;
% nsteps = 30;
% nsteps = 40;

%主循环
for n = 1:nsteps
    T = T+1;

    %计算D
    for j = 2:JE
        for i = 2:IE
            dz(i,j) = dz(i,j) + ...
                0.5 * ( hy(i,j) - hy(i-1,j) - hx(i,j) + hx(i,j-1) );
        end;
    end;

    % 加入高斯脉冲
    pulse = exp( -0.5 * (( t0-T ) / spread ) ^ 2  ) ;
    dz(ic,jc) = pulse;

    %计算E
    for j = 2:JE
        for i = 2:IE
            ez(i,j) = ga(i,j) * dz(i,j);
        end;
    end;

    %计算Hx
    for j = 1:JE-1
        for i = 1:IE-1
            hx(i,j) = hx(i,j) + 0.5*(ez(i,j) - ez(i,j+1));
        end;
    end;

    %计算Hy
    for j = 1:JE-1
        for i = 1:IE -1
            hy(i,j) = hy(i,j)  + 0.5 * (ez(i+1,j) - ..

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

不要迷恋我你会后悔的

离线夜的诱惑

人生若只如初见...

UID ：91870

注册：2012-04-11
登录：2012-06-03
发帖：46
等级：仿真新人

1楼发表于: 2012-05-06 21:57:04

clear all
%设置网格多少
IE=40;
JE=40;            % x,y,z方向的网格多少
KE=40;
%设置脉冲函数中心

ic=JE/2;
jc=JE/2;
kc=KE/2;
%初始化设置
ddx=0.01;         %空间网格大小
dt=ddx/6e8;       %时域网格大小
epsz=8.8e-12;     %介电常数
pi=3.14159;       %常数pi

%初始化系统变量
for k=1:KE
  for j=1:JE
     for i=1:IE
ex(i,j,k)=0.0;      %电场强度
ey(i,j,k)=0.0;
ez(i,j,k)=0.0;
dx(i,j,k)=0.0;      %电场通量D
dy(i,j,k)=0.0;
dz(i,j,k)=0.0;
hx(i,j,k)=0.0;      %x,y,z方向磁场强度
hy(i,j,k)=0.0;
hz(i,j,k)=0.0;
gax(i,j,k)=1.0;     %电场强度与D之间的关系
gay(i,j,k)=1.0;
gaz(i,j,k)=1.0;
     end;
  end;
end;

%指定偶极子
for k=11:30
   gaz(ic,jc,k)=0;
end;
   gaz(ic,jc,kc)=0;

%高斯脉冲变量设置
t0=20.0;
spread=6.0;
T=0;
nsteps=1;

while nsteps>0
  fprintf('nesteps-->');
  nsteps=input('Please input the number:');   %输入nsteps必须为整数
  fprintf('%d\n',nsteps);

for n=1:nsteps
  T=T+1;

%start the main fdtd loop
%calculate the Dx field
for k=2:KE
    for j=2:JE
        for i=1:IE
dx(i,j,k)=dx(i,j,k)+0.5*(hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1));
        end;
    end;
end;

%calculate the Dy field
for k=2:KE
    for j=1:JE
        for i=2:IE
dy(i,j,k)=dy(i,j,k)+0.5*(hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k));
        end;
    end;
end;

%calculate the Dz field
for k=1:KE
    for j=2:JE
        for i=2:IE
dz(i,j,k)=dz(i,j,k)+0.5*(hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k));
        end;
    end;
end;

%source
pulse=sin(2*pi*1500*1e6*T);
dz(ic,jc,kc)=pulse;
end;

%calculate the E from D field
for k=1:KE-1
    for j=1:JE-1
        for i=1:IE-1
ex(i,j,k)=gax(i,j,k)*dx(i,j,k);
ey(i,j,k)=gay(i,j,k)*dy(i,j,k);
ez(i,j,k)=gaz(i,j,k)*dz(i,j,k);
        end;
    end;
end;

%calculate the Hx field
for k=1:KE-1
    for j=1:JE-1
        for i=1:IE
hx(i,j,k)=hx(i,j,k)+0.5*(ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k));
        end;
    end;
end;

%calculate the Hy field
for k=1:KE-1
    for j=1:JE
        for i=1:IE-1
hy(i,j,k)=hy(i,j,k)+0.5*(ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k));
        end;
    end;
end;

%calculate the Hz field
for k=1:KE
    for j=1:JE-1
        for i=1:IE-1
hz(i,j,k)=hz(i,j,k)+0.5*(ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k));
        end;
    end;
end;

%end of the  main fdtd loop

%display
figure(1);

这里参照上面的Ez的图形怎么画  知道的指导下

end;

不要迷恋我你会后悔的

发帖回复