729阅读
0回复

[求助]恳请各位不吝赐教二维fdtd的程序解读

UID ：87816

注册：2012-01-10
登录：2022-01-12
发帖：56
等级：仿真一级

0楼发表于: 2012-05-03 18:37:34

本人刚刚接触fdtd，正在读这个程序，有些不明白的地方请大家指教
% 1.初始化
T=500; % 迭代次数
Nx=150;
Ny=150;
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*9*pi*10^9);
epsilon=1;
sigma=0;

ic=75; % 源的X位置
jc=75; % 源的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.333*(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;

      pulse=sin(2*pi*f*t*dt);  % 正弦波源
  %pulse=exp(-4*pi*(t*dt-t0)^2/(2/f)^2);%高斯脉冲情况

     ez(ic,jc)=ez(ic,jc)+pulse;  %软源，就是为了防止反射做的一个处理

    for i=2:Nx;  % 为了使每个电场周围都有磁场进行数组下标处理
        for j=2:Ny;
     ez(i,j)=gi3(i)*gj3(j)*ez(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1)); %
        end;
    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)=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;             ..

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

发帖回复

[求助]恳请各位不吝赐教 二维fdtd的程序解读

[求助]恳请各位不吝赐教二维fdtd的程序解读