[讨论]关于FDTD二维编程

%%我也来一个，是二维的 x?G"58  
clear;clc; OGnuBK
 

6"c(5#H  
%%%%%%%%%%%%%%% 1.初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |#O>DdKHT
 
rn-CQ2{?  
T=200;                                                          % 迭代次数 yT C+5_7  
Nx=100;  -=4:qQEw  
Ny=100; EUuk%<q7C(  
npml=8;                                                      %PML的网格数量 ,D{D QJ(B  
_VLA2#V>  
c0=3*10^8; PCiwQ4~  
f=1.5*10^(9);
AKUmh  
lambda=c0/f; J@(69&  
-pm%F8{T]  
wl=10; iq`caoi  
dx=lambda/wl;                                            %取lamda/10的精度 qL!pDZk  
dy=lambda/wl; 5wm(gF_t  
pi=3.14159; 6tBe,'*  
.QM>^(o$Z  
dt=dx/(2*c0);                                  %由于是二维情况同，取t=dx/2c0 }m.45n/  
t0=20*dt;                                      %为高斯脉冲的仿真作准备 sK|+&BC
 
"l-R|>6~  
epsz=1/(4*pi*9*10^9); OP\m~1  
epsilon=1; 7^]KQ2fF 8  
sigma=0; #HmZe98[%  
N9BfjT}  
ic=Nx/4;                                                        % 源的X位置 U~Rs?JmTdD  
jc=Ny/4;                                                        % 源的Y位置 p2(Z(V7*  
78t:ge eX  
for i=1:Nx+1; cFloaCz  
    for j=1:Ny+1; 9<1dps=c  
        dz(i,j)=0;                                            %z方向电位移 !,JV<(7k  
        ez(i,j)=0;                                              % z方向电场 *L#\#nh7  
        hx(i,j)=0;                                              % x方向磁场 8@d@T V!n&  
        hy(i,j)=0;                                              % y方向磁场 yey]#M[y  
        ihx(i,j)=0; ->a|  
        ihy(i,j)=0; ")eY{C  
        iz(i,j)=0;                                      % z方向求和中间参量 qPgny/(  
    end; (cm8x  
end; 9wP,Z"  
d&t,^Hj  
for i=2:Nx;  f2Z(hYH~  
    for j=2:Ny; 9%^O-8!  
        ga(i,j)=1; ?4YLt|sn  
    end; zEs:OOM  
end; _+}#  
D[p`1$E-1v  
%%%%%%%%%%%%%%%%%%%%PML参数的设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% E8b:MY  
for i=1:Nx; J42/S [Rt  
    gi2(i)=1; Ew*SA  
    gi3(i)=1; }o@Dsx5  
    fi1(i)=0; `E8D5'tt  
    fi2(i)=1; )i$KrN6  
    fi3(i)=1; BKd?%V8:Q  
end -5#cfi4^*  
r5'bt"K\>  
for j=1:Ny; tk!5"`9N  
    gj2(j)=1; t`E e/L%  
    gj3(j)=1; N<T@GQwkS  
    fj1(j)=0; \JIyJ8FleC  
    fj2(j)=1; kJ.0|l0  
    fj3(j)=1; DN;3VT.-  
end 6=g! Hs{  
^W,~  
b VcA#7 uA  
for i=1:npml+1;                                          %设置PML层中的参数 )H[Pz.'ah0  
    xnum=npml+1-i; ..UA*#%1  
    xn=0.33*(xnum/npml)^3;          %这里的0.333并不是严格的计算，而是经验值 dc,qQM  
    gi2(i)=1.0/(1+xn); E{{Kzr2$  
    gi2(Nx-1-i)=1/(1+xn); D7v_<  
    gi3(i)=(1-xn)/(1+xn); U,7}VdO  
    gi3(Nx-1-i)=(1-xn)/(1+xn); Mk!bmFZOZ  
o4nDjFhh  
    xn=0.25*((xnum-0.5)/npml)^3;              %这里的0.25和0.333也是一个道理 5odXT *n  
    fi1(i)=xn; 4XIc|a Aa  
    fi1(Nx-2-i)=xn; G]O5irsV  
    fi2(i)=1.0/(1+xn); KLgg([  
    fi2(Nx-2-i)=1/(1+xn); ;2l|0:  
    fi3(i)=(1-xn)/(1+xn); WGHf?G/s  
    fi3(Nx-2-i)=(1-xn)/(1+xn); ;={3H_{3  
end F $1f8U8  
lYhC2f m_  
for i=1:npml+1; )w,<XJhg`  
    xnum=npml+1-i; YpEH(tq  
    xn=0.33*(xnum/npml)^3; /^=8?wK  
    gj2(i)=1.0/(1+xn); ^JAp#?N^9  
    gj2(Ny-1-i)=1/(1+xn); |X;|=
.  
    gj3(i)=(1-xn)/(1+xn); 3_ko=& B$  
    gj3(Ny-1-i)=(1-xn)/(1+xn); 8\HzFB  
@IV,sze  
    xn=0.25*((xnum-0.5)/npml)^3; D4$2'h  
    fj1(i)=xn; % !Ih=DZ  
    fj1(Ny-2-i)=xn; =TJ9Gr/R&:  
Q+ZZwqyxD  
    fj2(i)=1.0/(1+xn); KRb'kW  
    fj2(Ny-2-i)=1/(1+xn); #O^%u,mJj  
    fj3(i)=(1-xn)/(1+xn); $) m$c5!  
    fj3(Ny-2-i)=(1-xn)/(1+xn); eD!mR3Ai@D  
end 8~ u/gM  
#M@~8dAH}M  
%%%%%%%%%%%%%%%%%%%2.迭代求解电场和磁场%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% z8cefD9F  
Ix+eP|8F  
for t=1:T;
1C(sBU"  
    for i=2:Nx;                  % 为了使每个电场周围都有磁场进行数组下标处理 B0D  
        for j=2:Ny; w$"^)EG,7  
            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)); #|=Q5"wU  
        end; z['2  
    end;                                                      % 电场循环结束 VE?Aa  
    4|;Ys-Q  
    pulse=sin(2*pi*f*t*dt);                                      % 正弦波源 in`|.#  
    %pulse=exp(-4*pi*(t*dt-t0)^2/(2/f)^2);                      高斯脉冲情况 1H
 \  
    dz(ic,jc)=dz(ic,jc)+pulse;            %软源，就是为了防止反射做的一个处理 Ae
jM\#>  
    ;2[),k  
    for i=1:Nx;                  % 为了使每个电场周围都有磁场进行数组下标处理 L5r02VzbD
 
        for j=1:Ny; u9N/9  
            ez(i,j)=ga(i,j)* dz(i,j);          %反映煤质的情况都是放到这里的 a] 7nK+N  
        end; BJKv9x1jK  
    end;                                                  % 电荷密度循环结束 c/E'GG%Q%  
    HYl~)O>  
    for j=1:Ny; Y-lTPR<Eq  
        ez(1,j)=0; &cwN&XBY  
        ez(Nx,j)=0; IYv.~IQO  
    end CV)K=Br5&_  
^G4@cR.An  
    for i=1:Nx; @D>qo=KPM  
        ez(i,1)=0; ojc.ykP$  
        ez(i,Ny)=0; /h8100  
    end % _nmv  
934@Z(aUH  
    for i=1:Nx;                  % 为了使每个磁场周围都有电场进行数组下标处理 R]3j6
\  
        for j=1:Ny-1;
N9 h|_ax  
            curl_e=ez(i,j)-ez(i,j+1); P=l 7m*m  
            ihx(i,j)=ihx(i,j)+fi1(i)*curl_e; `R!Q(rePx  
            hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j)); <Yg6=e  
        end; *Sps^Wl  
    end;                                                    % 磁场HX循环结束 T"1=/r$Ft  
G|u)eW  
    for i=1:Nx-1;                % 为了使每个磁场周围都有电场进行数组下标处理 $$i Gs6az  
        for j=1:Ny; 9 !$&1|,*  
            curl_e=ez(i+1,j)-ez(i,j); s"R5'W\U  
            ihy(i,j)=ihy(i,j)+fj1(j)*curl_e; %<fs \J^k  
            hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j)); [}yPy))A  
        end; MV]`[^xQ5  
    end;                                                    % 磁场HY循环结束 `\GRY @cg  
    mesh(ez); 6Ta+f3V  
    drawnow; {.ypZ8JU  
end;

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

关于三维的，用PML处理起来确实麻烦！ gu:vf/  
  特别是含有两三种介质的，大家有没有比较简单的程序 拿出来共享共享 ;"Aj80  
一起学习.

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

  我也在学

楼主好，使用一次

好东西，顶哦

继续顶，好心楼主

同关注这个问题

把主题地板那个，好好