[讨论]关于FDTD二维编程

%%我也来一个，是二维的 vPVA^UPNV  
clear;clc; ?89ZnH2/  
vYYLn9}5  
%%%%%%%%%%%%%%% 1.初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% L5"|RI}  
\Oku<5  
T=200;                                                          % 迭代次数 &HPzm6.3  
Nx=100;  6l{=[\.Xa  
Ny=100; Ik92='Z  
npml=8;                                                      %PML的网格数量 j3!]wolY  
,O9rL :?  
c0=3*10^8; AYC22(  
f=1.5*10^(9); =,*4:
TU  
lambda=c0/f; ,| \62B`  
?pT\Ft V  
wl=10; ^'g1? F$_  
dx=lambda/wl;                                            %取lamda/10的精度 R*C  
dy=lambda/wl; ?~mw  
pi=3.14159; [vIHYp  
tRqg')y  
dt=dx/(2*c0);                                  %由于是二维情况同，取t=dx/2c0 J!%c
HqR  
t0=20*dt;                                      %为高斯脉冲的仿真作准备 <lx~/3
<m  
Rwy<#9R[x  
epsz=1/(4*pi*9*10^9); w Wx,}=  
epsilon=1;  IX|2yu4  
sigma=0; `OY_v=}  
7[V6@K!Al[  
ic=Nx/4;                                                        % 源的X位置 O_PKS$sz{  
jc=Ny/4;                                                        % 源的Y位置 F-=W7 D:[c  
dMA"% R  
for i=1:Nx+1; "hwG"3n1  
    for j=1:Ny+1; VW7 ?{EL7  
        dz(i,j)=0;                                            %z方向电位移 ;'o:1{Y  
        ez(i,j)=0;                                              % z方向电场 c'9-SY1'~  
        hx(i,j)=0;                                              % x方向磁场 /r4QDwu  
        hy(i,j)=0;                                              % y方向磁场 1I#S?RSb  
        ihx(i,j)=0; (z[|\6O  
        ihy(i,j)=0; kUl:Yj=&  
        iz(i,j)=0;                                      % z方向求和中间参量 Qi_De '@  
    end; 4_`ss+gk  
end; B:YUb{CJ  
nUZ+N)*  
for i=2:Nx;  y6[^
I'kz  
    for j=2:Ny; V9v80e {n4  
        ga(i,j)=1; t/6t{*-w  
    end; zUw9  
end; 94=aVM\>>  
y.zS?vv2g  
%%%%%%%%%%%%%%%%%%%%PML参数的设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% {y{&tzZ  
for i=1:Nx; $X.X_  
    gi2(i)=1; 0FDfB;  
    gi3(i)=1; Unl6?_  
    fi1(i)=0; Q%xvS,oI  
    fi2(i)=1; LnN6{z{M  
    fi3(i)=1; Ib)>M`J
 
end QfI)+pf  
]}g;q*!J  
for j=1:Ny; G3
j&8[  
    gj2(j)=1; fbbk;Rq.'3  
    gj3(j)=1; UuV<#N)  
    fj1(j)=0; [l+1zt0w0  
    fj2(j)=1; sK#)wjj\^  
    fj3(j)=1; lP@/x+6tg  
end );LkEXC_'  
uT5sLpA|6  
))M; .b.D  
for i=1:npml+1;                                          %设置PML层中的参数 {r9fKA  
    xnum=npml+1-i; W_zv"c  
    xn=0.33*(xnum/npml)^3;          %这里的0.333并不是严格的计算，而是经验值 _T\~AwVc<  
    gi2(i)=1.0/(1+xn); 3^>D |  
    gi2(Nx-1-i)=1/(1+xn); -U6" Ce  
    gi3(i)=(1-xn)/(1+xn); 4{d`-reHg  
    gi3(Nx-1-i)=(1-xn)/(1+xn); RH,1U3?  
W$x'+t5H  
    xn=0.25*((xnum-0.5)/npml)^3;              %这里的0.25和0.333也是一个道理 ~Wq[H  
    fi1(i)=xn; 9zj^\-FA_l  
    fi1(Nx-2-i)=xn; @:'swO/\<  
    fi2(i)=1.0/(1+xn); U]!
D=+  
    fi2(Nx-2-i)=1/(1+xn); KrQ8//Ih  
    fi3(i)=(1-xn)/(1+xn); oG! S(95  
    fi3(Nx-2-i)=(1-xn)/(1+xn); {+V1>6  
end E(i<3U"4h[  
N'L3Oa\%  
for i=1:npml+1; e1e2W
k  
    xnum=npml+1-i; {DPobyvwFk  
    xn=0.33*(xnum/npml)^3; J%?
'Q{
  
    gj2(i)=1.0/(1+xn); D`;Q?fC  
    gj2(Ny-1-i)=1/(1+xn); V 0<>Xo%  
    gj3(i)=(1-xn)/(1+xn); *$ZLu jy7  
    gj3(Ny-1-i)=(1-xn)/(1+xn); ~E2xIhV  
#,dE)  
    xn=0.25*((xnum-0.5)/npml)^3; 1D /{Y  
    fj1(i)=xn;
fP41B  
    fj1(Ny-2-i)=xn; O -a`A.  
jk0Ja@8PK  
    fj2(i)=1.0/(1+xn); JW!SrM xF  
    fj2(Ny-2-i)=1/(1+xn); Qrt[MJ+#  
    fj3(i)=(1-xn)/(1+xn); aqTMOWyeu  
    fj3(Ny-2-i)=(1-xn)/(1+xn); 5feCA ,v7  
end DrD68$,QN  
-[kbHrl&  
%%%%%%%%%%%%%%%%%%%2.迭代求解电场和磁场%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LoHL}1BG-  
&n% 3rC5{  
for t=1:T; tHhA_  
    for i=2:Nx;                  % 为了使每个电场周围都有磁场进行数组下标处理
r5!I|E  
        for j=2:Ny; [.}-nAN  
            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)); =sG9]a<I  
        end; %N!2 _uk5  
    end;                                                      % 电场循环结束 L#'B-G4&y  
    NX?J  
    pulse=sin(2*pi*f*t*dt);                                      % 正弦波源 QI@
!QU$K&  
    %pulse=exp(-4*pi*(t*dt-t0)^2/(2/f)^2);                      高斯脉冲情况 8$fiq}a  
    dz(ic,jc)=dz(ic,jc)+pulse;            %软源，就是为了防止反射做的一个处理 ,!"\L~6  
    P)?)H]J"  
    for i=1:Nx;                  % 为了使每个电场周围都有磁场进行数组下标处理 $]/a/!d  
        for j=1:Ny; *KP 60T  
            ez(i,j)=ga(i,j)* dz(i,j);          %反映煤质的情况都是放到这里的 Xa._  
        end; e[t+pnRh  
    end;                                                  % 电荷密度循环结束 +zpmy3Q  
    o=`FGowF
 
    for j=1:Ny; sjLI^#a  
        ez(1,j)=0; Hu .e@7  
        ez(Nx,j)=0; 1mw<$'pm0  
    end H26'8e  
&e\A v.n@-  
    for i=1:Nx; $ctY#:;pV{  
        ez(i,1)=0; {1^9*  
        ez(i,Ny)=0; >$]SYF29  
    end Ju$vuE
O  
#D"fCVIS  
    for i=1:Nx;                  % 为了使每个磁场周围都有电场进行数组下标处理 W&[}-E8<Y  
        for j=1:Ny-1; L0{ehpvM  
            curl_e=ez(i,j)-ez(i,j+1); m:77pE&o  
            ihx(i,j)=ihx(i,j)+fi1(i)*curl_e; pq`Bg`c  
            hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j)); 8=^o2&  
        end; MtAD&+3$  
    end;                                                    % 磁场HX循环结束 m/"\+Hv  
wL]7d3t  
    for i=1:Nx-1;                % 为了使每个磁场周围都有电场进行数组下标处理 7d0E9t;W  
        for j=1:Ny; Zy2@1-z6  
            curl_e=ez(i+1,j)-ez(i,j); oWVlHAPj  
            ihy(i,j)=ihy(i,j)+fj1(j)*curl_e; fu/v1Nhm  
            hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j)); u 89u#gCAC  
        end; Az2HlKF"L  
    end;                                                    % 磁场HY循环结束 *yGOmi

 
    mesh(ez); &Y
2mLPB  
    drawnow; y%CaaK=V3  
end;

谢谢以上贴出程序的朋友！

关于三维的，用PML处理起来确实麻烦！ %k#Q)zWJ  
  特别是含有两三种介质的，大家有没有比较简单的程序 拿出来共享共享 k3se<NL[  
一起学习.

有没有做手机电磁辐射的，交流一下

  我也在学

楼主好，使用一次

好东西，顶哦

继续顶，好心楼主

同关注这个问题

把主题地板那个，好好