1073阅读
5回复

[求助]自己写的FDTD，有点问题

UID ：75963

注册：2011-04-20
登录：2011-06-03
发帖：14
等级：仿真新人

0楼发表于: 2011-04-25 16:58:55

要求：二维FDTD ，自由空间TM波，激励源是时谐场
程序如下：
clear;clc;
%设置网格大小
IE=100;
JE=100;
c=3e8;%自由空间的光速
%设置激励源中心
ic=IE/2;
jc=JE/2;
%激励源变量设置
f=10e9;
omg=2*pi*f;
T=1/f;
lmd=c/f;
%初始化设置
ddx=lmd/20;%空间步长
dt=ddx/(2*c);%时间步长
mu=4.0*pi*1.0e-7;%自由空间磁导率
epsl=1.0/(c*c*mu);%自由空间介电常数
sigma=0;%自由空间电导率
%初始化系统变量
Ez=zeros(IE,JE);
Hx=Ez;
Hy=Ez;
%设置时间步
Nsteps=150;%时间步为90时可画图，时间步为150时报：Warning: Axis limits outside float precision, use ZBuffer or Painters instead. %Not rendering
k1=dt/(mu*ddx);
k2=dt/(epsl*ddx);
%主循环
for N=1:Nsteps
   t=N*dt;
   %加入激励源
   if t<T/2
   sgn=0.5*(1-cos(omg*t))*sin(omg*t);%引入窗函数
   else
   sgn=sin(omg*t);
   end;
   Ez(ic,jc)=Ez(ic,jc)+sgn;
   %计算Hx
   for j=1:JE-1
   for i=1:IE-1
   Hx(i,j)=Hx(i,j)-k1*(Ez(i,j+1)-Ez(i,j));
   end;
   end;
   %计算Hy
   for j=1:JE-1
   for i=1:IE-1
   Hy(i,j)=Hy(i,j)-k1*(Ez(i+1,j)-Ez(i,j));
   end;
   end;
%计算Ez
   for j=2:JE
   for i=2:IE
   Ez(i,j)=Ez(i,j)+...
   k2*(Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1));
   end;
   end;
   %画图
   figure(1)
   mesh(Ez);
   drawnow