一个二维的FDTD程序
- % 本程序实现2维TM波FDTD仿真
- % 此程序用PML设置吸收边界条件
- % FDTD_2D_kongqi_PML
- % 仅喊有Ez,Hx,Hy分量
- clear;clc;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% 1.初始化
- T=200; % 迭代次数
- IE=100; %
- JE=100;
- npml=8; %PML的网格数量
- c0=3*10^8; % 波速
- f=1.5*10^(9); % 频率
- lambda=c0/f; % 波长
- wl=10;
- dx=lambda/wl;
- dy=lambda/wl;
- pi=3.14159;
- dt=dx/(2*c0); % 时间间隔
- epsz=1/(4*pi*9*10^9); % 真空介电常数
- epsilon=1; % 相对介电常数
- sigma=0; % 电导率
- spread=6; % 脉冲宽度
- t0=20; % 脉冲高度
- ic=IE/2; % 源的X位置
- jc=JE/2; % 源的Y位置
- for i=1E+1;
- for j=1:JE+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=2E; %
- for j=2:JE;
- ga(i,j)=1;
- end;
- end;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %PML参数的设置
- for i=1E;
- gi2(i)=1;
- gi3(i)=1;
- fi1(i)=0;
- fi2(i)=1.0;
- fi3(i)=1.0;
- end
- for j=1:JE;
- 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;
- gi2(i)=1.0/(1+xn);
- gi2(IE-1-i)=1/(1+xn);
- gi3(i)=(1-xn)/(1+xn);
- gi3(IE-1-i)=(1-xn)/(1+xn);
- xn=0.25*((xnum-0.5)/npml)^3;
- fi1(i)=xn;
- fi1(IE-2-i)=xn;
- fi2(i)=1.0/(1+xn);
- fi2(IE-2-i)=1/(1+xn);
- fi3(i)=(1-xn)/(1+xn);
- fi3(IE-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(JE-1-i)=1/(1+xn);
- gj3(i)=(1-xn)/(1+xn);
- gj3(JE-1-i)=(1-xn)/(1+xn);
- xn=0.25*((xnum-0.5)/npml)^3;
- fj1(i)=xn;
- fj1(JE-2-i)=xn;
- fj2(i)=1.0/(1+xn);
- fj2(JE-2-i)=1/(1+xn);
- fj3(i)=(1-xn)/(1+xn);
- fj3(JE-2-i)=(1-xn)/(1+xn);
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% 2.迭代求解电场和磁场
- for t=1:T;
- for i=2E; % 为了使每个电场周围都有磁场进行数组下标处理
- for j=2:JE;
- 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); % 正弦波源
- dz(ic,jc)=dz(ic,jc)+pulse; % 软源
- for i=1E; % 为了使每个电场周围都有磁场进行数组下标处理
- for j=1:JE;
- ez(i,j)=ga(i,j)* dz(i,j); %反映煤质的情况都是放到这里的
- % iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ;
- end;
- end; % 电荷密度循环结束
- for j=1:JE;
- ez(1,j)=0;
- ez(IE,j)=0;
- end
- for i=1E;
- ez(i,1)=0;
- ez(i,JE)=0;
- end
- for i=1E; % 为了使每个磁场周围都有电场进行数组下标处理
- for j=1:JE-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=1E-1; % 为了使每个磁场周围都有电场进行数组下标处理
- for j=1:JE;
- 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循环结束
- end;
- end;
复制代码
[ 本帖最后由 风花雪月 于 2006-11-12 07:39 编辑 ] |