• love12121
    了解作者
  • matlab
    开发工具
  • 1KB
    文件大小
  • rar
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 18
    下载次数
  • 2009-11-24 17:00
    上传日期
Matlab写的二维带PML边界激励源是高斯脉冲的FDTD程序
Prog_2D_FDTD.rar
  • Prog_2D_FDTD.m
    3.4KB
内容介绍
%高斯脉冲入射条件下的FDTD模拟程序 clear;clc; % 1.初始化 T=2000; % 迭代次数 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; 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; end; % 磁场HX循环结束 for i=1:Nx-1; % 为了使每个磁场周围都有电场进行数组下标处理 for j=1:Ny; curl_e=ez(i+1,j)-ez(i,j); ihy(i,j)=ihy(i,j)+fj1(j)*curl_e; hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j)); end; end; % 磁场HY循环结束 %mesh(ez); pcolor(abs(ez)); %drawnow; shading flat; pause(0.002); end;
评论
    相关推荐