[讨论]关于FDTD二维编程

%%我也来一个，是二维的 Ee2P]4_d  
clear;clc; 2c 0;P #ol  
AX8~w(sv  
%%%%%%%%%%%%%%% 1.初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RI&O@?+U  
 4d5c]%  
T=200;                                                          % 迭代次数 FL^ _)`  
Nx=100;  @6>Q&GYqt  
Ny=100; E6TeZ%g  
npml=8;                                                      %PML的网格数量 j'z#V_S  
R+Hu?Dv&F  
c0=3*10^8; WJhTU@'  
f=1.5*10^(9); 7,+eG">0  
lambda=c0/f; <v|"eq}  
.iRKuBM/  
wl=10; &/]en|f"  
dx=lambda/wl;                                            %取lamda/10的精度 1fS&KO{a  
dy=lambda/wl; k773h`;  
pi=3.14159; cZNcplt8  
r6Yd"~
n  
dt=dx/(2*c0);                                  %由于是二维情况同，取t=dx/2c0 E=ijt3  
t0=20*dt;                                      %为高斯脉冲的仿真作准备 e=OHO,74z"  
+9b{Y^^~T  
epsz=1/(4*pi*9*10^9); a31
e.36g  
epsilon=1; )bOfs*S  
sigma=0; @/.# /  
DJ,LQj  
ic=Nx/4;                                                        % 源的X位置 uI7n{4W*x  
jc=Ny/4;                                                        % 源的Y位置 0]|`*f&p;  
='o3<}  
for i=1:Nx+1; ]u|v7}I4  
    for j=1:Ny+1; <J&S[`U!  
        dz(i,j)=0;                                            %z方向电位移 _M%S  
        ez(i,j)=0;                                              % z方向电场 sX%n`L  
        hx(i,j)=0;                                              % x方向磁场 Q:(mK* _  
        hy(i,j)=0;                                              % y方向磁场 pD%(Y^h?  
        ihx(i,j)=0;  n?EgC8b9  
        ihy(i,j)=0; O@T,!_Zf  
        iz(i,j)=0;                                      % z方向求和中间参量 ^[xcfTN  
    end; CW &z?Bra  
end; hTX[W%K  
Kq3c Kp4  
for i=2:Nx;  noxJr/A]  
    for j=2:Ny; 51q|-d  
        ga(i,j)=1; \Mg_Q$  
    end; <#C,66k  
end; 8@m$(I+  
H~+A6g]T  
%%%%%%%%%%%%%%%%%%%%PML参数的设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% UPuoIfuqI  
for i=1:Nx; v.-DXQ
q  
    gi2(i)=1; [P`e@$  
    gi3(i)=1; CLrX!JV>  
    fi1(i)=0; S\).0goOW  
    fi2(i)=1; 9;e!r DW,#  
    fi3(i)=1; d2`g,~d  
end -+rzc&h  
\sAaVdZJH(  
for j=1:Ny; 5)!g.8-!  
    gj2(j)=1; Ik0g(-d  
    gj3(j)=1; \FVfV`x  
    fj1(j)=0; $9,&BW_*  
    fj2(j)=1; A9MM^jV8  
    fj3(j)=1; )RJEOl1  
end \~ D(ww  
ZZT #V%Q=u  
H{*Dc_  
for i=1:npml+1;                                          %设置PML层中的参数 F}nwTras  
    xnum=npml+1-i; p@Y=6Bw  
    xn=0.33*(xnum/npml)^3;          %这里的0.333并不是严格的计算，而是经验值 y!#-[K:  
    gi2(i)=1.0/(1+xn); 4b2d(x)0X  
    gi2(Nx-1-i)=1/(1+xn); '>"{yi-  
    gi3(i)=(1-xn)/(1+xn); B}iEhWO6  
    gi3(Nx-1-i)=(1-xn)/(1+xn); }<04\t?  
1hG#  
    xn=0.25*((xnum-0.5)/npml)^3;              %这里的0.25和0.333也是一个道理 )!"fUz$  
    fi1(i)=xn; >NKJ@4Y  
    fi1(Nx-2-i)=xn; AoS7B:T;!  
    fi2(i)=1.0/(1+xn); rm[C{Pn  
    fi2(Nx-2-i)=1/(1+xn); )\akIA  
    fi3(i)=(1-xn)/(1+xn); U g"W6`  
    fi3(Nx-2-i)=(1-xn)/(1+xn); ,dIev<  
end HT;^
u"a~  
XIdh9)]^}  
for i=1:npml+1; .G O0xn
m  
    xnum=npml+1-i; 8q[; 0  
    xn=0.33*(xnum/npml)^3; Tlw'05\{J  
    gj2(i)=1.0/(1+xn); 7Z6=e6/\  
    gj2(Ny-1-i)=1/(1+xn); puC91  
    gj3(i)=(1-xn)/(1+xn); :>1nkm&Eg  
    gj3(Ny-1-i)=(1-xn)/(1+xn); yq%5h[M  
L^sjV/\oW  
    xn=0.25*((xnum-0.5)/npml)^3; *jlIV$r_  
    fj1(i)=xn; YMX9Z||  
    fj1(Ny-2-i)=xn; 4@PA+(kvS  
h[mT4e3c  
    fj2(i)=1.0/(1+xn); yH/A9L,Z  
    fj2(Ny-2-i)=1/(1+xn); .Kb3VNgwvm  
    fj3(i)=(1-xn)/(1+xn); :o'x?]  
    fj3(Ny-2-i)=(1-xn)/(1+xn); L'= \|r  
end &xnQLz:#  
.{V"Gn9!  
%%%%%%%%%%%%%%%%%%%2.迭代求解电场和磁场%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ~bLx2=-"  
3kn-tM  
for t=1:T; QKG3>lU  
    for i=2:Nx;                  % 为了使每个电场周围都有磁场进行数组下标处理 %j7b0pb  
        for j=2:Ny; \vA*dQ-  
            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)); <Y]LY_(  
        end; 3:+9H}Q  
    end;                                                      % 电场循环结束 3z8i0  
    1t
 R^  
    pulse=sin(2*pi*f*t*dt);                                      % 正弦波源 +hyOc|5  
    %pulse=exp(-4*pi*(t*dt-t0)^2/(2/f)^2);                      高斯脉冲情况 4ijtx)SA  
    dz(ic,jc)=dz(ic,jc)+pulse;            %软源，就是为了防止反射做的一个处理 YWFE*wQ!  
    Hb)FeGsd).  
    for i=1:Nx;                  % 为了使每个电场周围都有磁场进行数组下标处理 %fj5;}E.  
        for j=1:Ny; c7e,lgG-  
            ez(i,j)=ga(i,j)* dz(i,j);          %反映煤质的情况都是放到这里的 }Pcm'o_wT  
        end; Pxlc RF  
    end;                                                  % 电荷密度循环结束 #a
8kA"X  
    0pZvW  
    for j=1:Ny; PF%-fbh!~  
        ez(1,j)=0; :7v'[b  
        ez(Nx,j)=0; M-o'`e'  
    end IVR%H_uz  
y)G-6sZ/  
    for i=1:Nx; V@ >(xe7  
        ez(i,1)=0; n>Ff tVZNJ  
        ez(i,Ny)=0; y{ur'**l  
    end k#].nQG  
KF.{
r  
    for i=1:Nx;                  % 为了使每个磁场周围都有电场进行数组下标处理 8SGaS&  
        for j=1:Ny-1; :X0L6y)u  
            curl_e=ez(i,j)-ez(i,j+1); w*qj0:i5as  
            ihx(i,j)=ihx(i,j)+fi1(i)*curl_e; /I%z7f91O  
            hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j)); L+am-k:T~  
        end; @5>#<LV=E#  
    end;                                                    % 磁场HX循环结束 [X(4
( 1i  
nKx)R^]k  
    for i=1:Nx-1;                % 为了使每个磁场周围都有电场进行数组下标处理 b#VtPn]  
        for j=1:Ny; (>al-vZ6A  
            curl_e=ez(i+1,j)-ez(i,j); `b?uQ\#-M  
            ihy(i,j)=ihy(i,j)+fj1(j)*curl_e; UqP{Cyy{  
            hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j)); ^hIdmTf6  
        end; M8/a laoT  
    end;                                                    % 磁场HY循环结束 Zy}tZRG  
    mesh(ez); 7~7L5PRW  
    drawnow; *1)>He$qL  
end;

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

关于三维的，用PML处理起来确实麻烦！ cAWn*%  
  特别是含有两三种介质的，大家有没有比较简单的程序 拿出来共享共享 (@?mm  
一起学习.

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

  我也在学

楼主好，使用一次

好东西，顶哦

继续顶，好心楼主

同关注这个问题

把主题地板那个，好好