2d-fdtd.m

微波仿真论坛_高斯脉冲入射条件下的FDTD模拟程序

clear;clc;  
 
% 1.初始化 
 
T=200; % 迭代次数   
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;  
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/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);  
    drawnow;  
end;