FDTD_面波模拟_模拟平面波_fdtd_matlab_FDTD平面波matlab

clear;clc; %%e799bee5baa6e58685e5aeb931333262343133%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IE=100; JE=100;% 网格数量 L=400; % 迭代次数 f=1.6e9; % 频率 dx=.01; % 网格尺寸 dt=dx/2/3e8; % 时间间隔 ic=3;%平面波源位置 e0=8.8e-12; % 真空介电常数 m0=pi*4e-7; % 真空磁导率 epsilon=1; % 相对介电常数 sigma=0; % 电导率 ia=7; ib=IE-ia+1; ja=7; jb=JE-ja+1; for i=1:IE for j=1:JE ga(i,j)=1/(epsilon+sigma*dt/e0); % 求和参量 gb(i,j)=sigma*dt/e0; % 求和参量 end end ez_inc=zeros(1,IE); hx_inc=zeros(1,JE); ez_low1=0; ez_low2=0; ez_low3=0; ez_low4=0; dz=zeros(IE,JE); % z方向电荷密度 ez=zeros(IE,JE); % z方向电场 iz=zeros(IE,JE); % z方向电场求和参量 hx=zeros(IE,JE); % x方向磁场 hy=zeros(IE,JE); % y方向磁场 ihx=zeros(IE,JE); % x方向磁场参量 ihy=zeros(IE,JE); % y方向磁场参量 %%%%%%%%%%%%%%%%%%%%%%%%%迭带求解电磁场%%%%%%%%%%%%%%%%% for T=1:L for j=2:JE ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j)); end ez_inc(1)=ez_low2; ez_low2=ez_low1; ez_low1=ez_inc(2); ez_inc(JE)=ez_low3; ez_low3=ez_low4; ez_low4=ez_inc(JE-1); %%%%%%%%%%%%%% 电荷密度dz%%%%%%%%%%%%%%%%%%%%%%%%% for i=2:IE-1 for j=2:JE dz(i,j)=dz(i,j)+0.5*( hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1)); end end source=sin(2*pi*f*dt*T); ez_inc(ic)=source; for i=ia:ib dz(i,ja)=dz(i,ja)+0.5*hx_inc(ja-1); dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb); end %%%%%%%%%%%%%% 电场ez%%%%%%%%%%%%%%%%%%%%%%%% for i=1:IE for j=1:JE ez(i,j)=ga(i,j)*( dz(i,j)-iz(i,j) ); iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ; end end %%%%%%%%%%%%%%%%%%%%边界电场%%%%%%%%%%%%%%%%%%%%%%% for i=1:IE ez(1,j)=0; ez(IE,j)=0; end for i=1:JE ez(i,1)=0; ez(i,JE)=0; end for j=1:JE-1 hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1)); end %%%%%%%%%%%%%%%%%%%%%%%%磁场Hx%%%%%%%%%%%%%%%%%%%%%% for i=1:IE for j=1:JE-1 hx(i,j)=hx(i,j)+0.5*(ez(i,j)-ez(i,j+1)); end end for i=ia:ib hx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja); hx(i,jb)=hx(i,jb)-0.5*ez_inc(jb); end %%%%%%%%%%%%%%%%%磁场Hy%%%%%%%%%%%%%%%%%%% for i=1:IE-1 for j=1:JE hy(i,j)=hy(i,j)+0.5*(ez(i+1,j)-ez(i,j)); end end for j=ja:jb; hy(ia-1,j)=hy(ia-1,j)-0.5*ez_inc(j); hy(ib,j)=hy(ib,j)+0.5*ez_inc(j); end %%%%%%%%%%%%%%%%%%绘图%%%%%%%%%%%%%%%%%% surf(ez) axis([0 100 0 100 -2 2]) shading interp; %view(0,90); drawnow; end